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the propagation of shock waves through regions of 
non-uniform area or flow 
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SUMMARY 

This paper refers to the work of Moeckel (1952) on the 
interaction of an oblique shock wave with a shear layer in steady 
supersonic flow and the work of Chester (1955) and Chisnell (1957) 
on the propagation of a shock wave down a non-uniform tube. 
It is shown that their basic results can be obtained by the 
application of the following simple rule. ‘The relevant equations 
of motion are first written in characteristic form. ‘Then the rule 
is to apply the differential relation which must be satisfied by the 
flow quantities along a characteristic to the flow quantities just 
behind the shock wave. ‘Together with the shock relations this 
rule determines the motion of the shock wave. The accuracy of the 
results for a wide range of problems and for all shock strengths is 
truly surprising. 

The results are exactly the same as were found by the authors 
cited above. The derivation given here is simpler to perform 
(although the original methods were by no means involved) and 
of somewhat wider application, but the main reason for presenting 
this discussion is to try to throw further light on these remarkable 
results. 

In discussing the underlying reasons for this rule, it is 
convenient to use the propagation in a non-uniform tube as a 
typical example, but applications to a number of problems are 
given later. A list of some of these appears at the beginning of 
the introductory section. 


1. INTRODUCTION 


There are various problems arising in the different branches of fluid 


dynamics which have as a main feature the interaction of a shock wave 


with a non-uniform region of some sort. Some examples are: 


(i) the motion of a shock wave down a non-uniform channel or tube; 
(ii) the propagation of a shock normally through a plane distribution 
of density, pressure, etc. ; 
(iii) converging cylindrical and spherical shocks (including non-uniform 
states ahead of the shock such as occur in magnetohydrodynamics) ; 
F.M. Y 
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(iv) the interaction of an oblique shock with a shear layer in steady 
supersonic flow ; 

(v) the shape and strength of the shock inside the entry of an axi- 
symmetrical supersonic duct; 

(vi) the propagation of a bore in water of non-uniform depth; 

(vii) the motion of kinematic shock waves in traffic flow and flood waves 

when the flow conditions downstream are non-uniform. 
We shall ultimately consider all of these problems, but in the general 
discussion we may refer to (i) as a typical example. It should be stated at 
once that only the one-dimensional (hydraulic) formulation is intended 
in (i), i.e. all quantities are averaged over the cross-sectional area and become 
functions of the time and distance down the tube only. 

Of course, even in uniform media, the problem of shock propagation 
cannot be solved analytically in general and numerical methods have to 
be used. But to study some of the main effects of the interaction we can 
consider cases in which changes in the speed and strength of the shock 
are caused entirely by the non-uniform region. ‘Thus, for the propagation 
down a non-uniform tube we suppose that to the left (~ < 0) of some 
section x = 0, the cross-sectional area A(x) is constant and the shock 
wave initially moves in this constant part of the tube travelling with 
uniform speed. ‘Taking the usual model, we can suppose that the shock 
wave is caused by a piston moving with constant speed; the piston is 
always taken far enough to the left for the reflections from it not to come 
back and interfere with the basic interaction. ‘The analogous suppositions 
for the other problems listed above are obvious except perhaps in the case 
of (iil) and (v). But the converging shocks of (111) can be treated as special 
cases of the non-uniform tube. The cylindrical shock is obtained by 
choosing a wedge-shaped channel with A(x) < (x)—.x) and the spherical 
shock by taking a cone-shaped tube with A(x) x (xy—.x)*. In_ this 
formulation we can still suppose there is a uniform section with A constant 
before the varying part. However, Guderley’s similarity solutions for 
converging shocks (Guderley 1942) will be discussed below and they do 
not correspond exactly to these initial conditions. But converging 
cylindrical shocks are known to be insensitive to the details of the initial 
conditions. For example, Payne’s numerical calculations (Payne 1957) 
agree closely with Guderley’s solution although the initial conditions are 
quite different (in fact, they are close to the ones proposed above). ‘The 
same thing will be found here, with similar results for problem (v). 

Apart from the special features of (iii) and (v) the initial conditions 
adopted above will not be exactly fulfilled in practice in the other cases; 
even without the non-uniformities in the medium, the shock may be varying 
due to changes in the flow conditions behind it (such as variations in the 
speed of the effective piston in the model mentioned above). But in many 
cases, the scale of this variation will be much greater than that of the 
interaction and so will have a negligible effect on the ‘local’ behaviour in 


the non-uniform region. 
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‘The considerations here arise from the work of Moeckel (1952), Chester 
(1954), and Chisnell (1955, 1957). Moeckel studied problem (iv) by the 
following method. ‘lhe shear layer is represented by a set of parallel 
surfaces at which small discontinuities in velocity, pressure, density, etc., 
occur. ‘lhe interaction of an oblique shock with any of these can be solved 
to find the resulting change in the strength of the shock in terms of the 
change in the incident flow. ‘The successive interactions are not inde- 
pendent, but by ignoring the dependence between them comparatively 
simple laws are obtained tor the change in strength of the shock as it 
progresses through the shear layer. With this basic result Moeckel goes 
on to discuss refinements. Chisnell uses essentially the same methods to 
study problems (i) and (11); on the whele, it is more convenient to consider 
these as far as basic ideas are concerned. In the first of them the elementary 
interaction is between the shock and a small change 6/4 in cross-sectional 
area; the resulting change 5.V in the Mach number ™ of the shock wave 
is given by the formula 

dA 2MbM \ 
A (M®-1)K(M)’ (1) 
where K() is a slowly varying function which starts at 0-5 for weak shocks, 
M = 1, and tends to 0-394 (for y = 1-4) as M— «. Chisnell represents 
the tube by a series of such small changes and for a first approximation 





neglects the interference between successive interactions. ‘This amounts to 
integrating (1) to get _Wasafunctionof A. Inthelimitas WM -> 1, this gives 


M-1« A-, (2) 
which is the law of geometrical acoustics; in the limit as M—-> =m, it gives 
M x A-'K x, (3) 


where K,, = 0-394 for y = 1-4, as noted above. Now (2) is well established 
for weak shocks* and Chisnell checked (3) by applying his formula to 
converging cylindrical and spherical shocks and comparing with Guderley’s 
exact solutions for these cases. ‘The cylindrical shock corresponds to a 
wedge-shaped tube with 4 proportional to the distance from the axis and 
the spherical one corresponds to a conical tube with 4 proportional to the 
square of the distance from the centre. ‘Thus, according to (3), the exponents 
in the formula for MW asa function of distance are }K , and K_, respectively. 
The comparison with the Guderley values (as calculated by Butler (1954)) 
is shown in the following table: 








; 
. | 
y Cylindrical shock | Spherical shock 
| oe ! } | 

| Chisnell | Guderley | Chisnell | Guderley | 

| | | s 
|} 6/5 0163112 | 0-161220 | 0:326223 | 0-320752 | 
7/5 0197070 | 0:197294 | 0394141 | 0:394364 | 
5/3 0-225425 0-226054 0-452108 | 0-452692 | 
| 


| = 4 | 


| 1 i 





In the present circumstances this is true; there are important differences between 
weak shocks and acoustic pulses in general. 


v2 
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Comment on these is unnecessary! Payne (1957) determined the motion 
of converging cylindrical shocks numerically for a whole range of shock 
strengths and found that in all cases Chisnell’s law was extremely accurate. 
Chisnell went on to find a higher approximation by including re-reflected 
waves and the corrections turned out to be of the expected small magnitude. 
In his approach the high accuracy of the first approximations appears to 
some extent as a coincidence. 

Chester’s contribution was on different lines. He considered problem (i) 
for the case in which the cross-sectional area A(x) remains close to some 
mean value. He developed a theory for the small perturbations in the flow 
behind the shock and solved the linearized equations exactly. The formula 
(1) was first obtained by Chester in this way. Although Chester’s work 
is more restricted than Chisnell’s, it adds greatly to an understanding of 
these results; it is discussed in detail in the next section. 

Now, the purpose of the present paper is to point out that the basic 
formulae obtained by these authors can be derived in the following significant 
way. First the appropriate equations of motion for the flow are written 
in characteristic form. For example, for a non-uniform tube we have 

dp +pa du+ eee 0 (4) 

u+aA 

ona positive characteristic dx/dt = u +a, where p, p, a, u denote the pressure, 
density, sound speed, particle velocity, respectively. Then, the character- 
istic relation is applied (quite illogically it may seem) to the flow quantities 
at the shock wave. But these quantities are all known in terms of the shock 
strength from the Rankine—-Hugoniot shock relations. ‘Thus on substitution 
in the characteristic equation, an equation for the variation of the shock 
strength is obtained. For example, in (4) we substitute the expressions 
given by the shock relations for p, p, a, u in terms of the Mach number of 
the shock. ‘This gives a first-order equation for M as a function of A which 
can be integrated immediately. It is exactly Chisnell’s formula, the 
differential equation for M(A) being identical with (1). The constan 
of integration is fixed from the initial strength of the shock in the straight 
portion of the tube. When the area A is constant the variations in shock 
strength caused by non-uniform conditions ahead of the shock come in 
only through the shock relations. 

Even for quite complicated equations it is usually a simple matter to 
write down the characteristic equation corresponding to (4). In addition 
only the shock relations are required—and these are needed in any approach. 

Once this derivation has been noticed it is easy to see why it agrees 
with the previous methods. Full details are given in the subsequent sections. 
New questions also arise. For instance, this new approach is found to have 
a close connection with Butler’s method for calculating the Guderley 
solutions for converging shocks. Again the method is related in some way 
with ‘shock-expansion theory’. These various topics and additional 
applications are discussed below; the section headings are self-explanatory. 
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In concluding these general remarks it should be said that the discussions 
in this paper still fall short of a full understanding of all the questions 
involved; it is still not completely clear to what extent the unexpected 
accuracy of the results in some cases should be ascribed to coincidence. 

2. DETAILED DISCUSSION OF THE METHOD 

As mentioned already, we shall consider first the typical example of the 
propagation of a shock down a tube of variable area A(x) containing gas 
originally at rest with constant pressure p) and density py. We suppose 
that A(x) takes a constant value 4, in x < 0 and that the shock is initially 





Shock 





Figure 1. The (x, t)-plane for the interaction of a shock wave with a non-uniform 
region. 


moving with constant velocity UU, in this region. ‘The flow quantities 
Pi» Pr» 4, U, behind the shock are determined in terms of U, by the shock 
conditions. Conversely, given u, (the velocity of a piston maintaining the 
motion for example) the shock velocity U, and the other quantities can 
be found. When the shock reaches x = 0, disturbances are propagated 
back into the uniform flow region as represented in figure 1, and the future 
motion of the shock is modified. This reflected disturbance propagates 
along negative characteristics labelled C_ in the figure; in addition, entropy 
changes are carried along the particle paths labelled P. From the physical 
point of view the positive characteristics (one labelled as Cis shown in 
figure 1) play a subsidiary role. We may note that in xv < 0 the C_ form 
a simple wave, 1.e. they are all straight and 2a/(y—1)+u = 2a,/(y—1)+u, 
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throughout. 


G 
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(We assume that if a reflected shock wave is formed it does 


so outside the region under consideration. ) 
‘The equations of motion are 


Ppt Up, t a(u, t 


Uy 





> a 


1 
+uu,+ -p, = 9, 
p 


Pr Tr up, —a*(p,+up,) =U, | 


where a® yp /p- 


‘lhe characteristic form (1.e 


. the form for which each 


equation contains derivatives in only one direction in the (.x, t)-planc) 


is 
p u, pauA 
t_ +p.) tpa|—— +2, ) + ———. = ©, | 
uta uta (uta)A 
(6) 
Pr of Pt 
—+p,})—-—a([—+p,] =. 
u Ou 
We can also write these as 
au dA dx 
dp t+ pa lu . = a oY on & — u - 
ikl ice ita A dt — 
tu dA dx 
lp — pa du vas ao =) gS sas u—a, (7 
= aa dt my 
' 1x 
dp—a*dp = on P: — = tl. 
dt 


The rule is to apply the tirst of these (valid along a C,) to the shock wave. 


From the shock conditions 











? 


{ 2 ee | 
p py as {‘— Ve - ———~}, 
iy + Vly +1) 
(yv+1)M* 
f diy ae a =. Pe) 
Pot —1)M242 °) 
Pa; l 
“=a, a” a) 
p+] M 
where the shock velocity | a,M. Substituting these expressions into 
the first of (7), we have 
2M dM dA 
sree ti p= 6): (9) 
(.V?—1) K(M) A 
where 
2 1-p2 ; 
K(M) 2| (1 —— )Qu+ tear) | ' (10) 
‘+1 op 
% (< u\? (y—1)M?+2 
‘ a ) 2yM?—-(y—1) 
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The law of propagation of the shock is given by the function M(4) which 
satisfies (). 
As M +1, » +1 and K -- 0-5; hence, the solution of (9) is 


M-1a A+ 


= (394 for y = 1-4; (11) 
hence, 
M x A-'8 a, 
‘he corresponding laws for the variations of p, p, w are then found from (8). 
We must now investigate why this rule works so well. For some reason 


pau dA 
i+aA 





dp +pa du + 


is almost zero along the shock. ‘This means that 


A PP de 
(F Pe) pa( 7 u,) utaA 


is very small at the shock. But, by making use of the first of equations (6) 
(which holds everywhere), we conclude that 


l ( 
(Z ~ aa) (p; + pau, ) (12) 





is very small at the shock. Now the obvious first thought as to why the 
rule works is that the characteristic C. follows along behind the shock 
very close to it (see figure 1), and so the relation valid along the C, is a 
good approximation along the shock. ‘This would appeal to a supposed 
smallness of the first factor in (12), but although (w+a—U)/U is zero 
for M = 1, it tends to [4/{2y(y—1)} —(y—1)]/(y +1) = 0-274 (for y = 1-4) 
as M -+ », which is at least a hundred times the error found in the results. 
In fact, it is the smallness of the second factor in (12) which leads to the high 
accuracy. 

In Chester’s small perturbation theory it will be shown that p, + pau, is 
zero at the shock; in fact p,+ pau, is zero everywhere and the above result 
is the correct answer in that theory. However, when finite changes in A 
are considered, as in Chisnell’s work on the cylindrical shock, it is not clear 
why this should be so. From the numerical calculations it is observed that 
p, + pau, is very small, but all the analytical approaches considered so far 
work on expansions assuming essentially that the first factor in (12) is small. 
Then the high accuracy of the lowest order results comes out as something 
of a coincidence. However, these approaches are very relevant to Butler’s 


calculations and some comments are made in the next section. 
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First we consider Chester’s small perturbation theory. It is clear at 
the outset that if p, p, a, uw remain close to the original values p,, p,, a), uy 
respectively, then in the linearized theory we have 
py Qj uy, dA 


dp + py @y du + =() 


u,+a,A, 
oneach C,. Since each C, starts in the region where p =. p, etc., it follows 
that 

py au, A(x) — A, 


u,+a, A, 





(p—Pi) +p. a (u-—u) + i (13) 
throughout the flow and in particular at the shock. On substituting the 
shock relations the result derived above is obtained. Also, differentiating (13) 
with respect to ¢, we see that p,+p,a@,u, = 0 everywhere. ‘The relation (13) 
must also apply across one of Chisnell’s interactions and this explains why 
the present method gives the same answer as Chisnell’s. 

It is interesting to look at Chester’s full solution. ‘The equations (5) 
may be linearized in the straightforward way* and the general solution 
is found to be 


p, aj uz A(x)—A 











) L + F(x—ju,+a,'t)+ G(x—{ 't), (14 
P=Pi w—ae A, + F(x— {uy +a,jt) + G(x—{u,—ayjt), (14) 
PF el ens ee ee 

ut — a? A, pi a 
‘ G(x—{u,—a,}t), (15) 
Pi 
pP— py aoe + H(x—u,t), (16) 


where F, G and H are arbitrary functions. ‘This solution shows very 
clearly how the disturbances are made up of the four contributions; first, 
terms directly from the area change, then disturbances propagating on 


v—{u, +a,}t = constant (C,), v—{u,—a,{t = constant (C_) 
and v—u,t = constant (P). 


The entropy changes are represented by the function // and it is interesting 
to note that p and w do not depend directly on them. ‘lhe shock conditions 
provide two relations between p, p and u; these boundary conditions serve 
to determine the functions G and H. ‘The other function F is determined 
by the flow in x < 0 and in the present case F = OF. With F = 0 it is 
observed from (14) and (15) that p,+p, a, u, = 0. 


* Actually Chester gives a more thorough treatment: he starts with the full three- 
dimensional equations of motion and carries through in detail an averaging process 
which leads to the hydraulic theory. 

+ Contributions to F arise if the piston motion in the straight portion of the tube 
is not uniform. 
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In general, p,+p,a,u, = —2(u,+a,)F(x— {uy +a,}t) so that p, +p, a, u; is 
constant along a positive characteristic C ,; the t-derivatives knock out the 
terms depending directly on A and this particular combination gets rid 
of the G. In the general problem when the disturbances are not small, 
we may expect that in any small region the small perturbation theory holds 
for variations about some local values. ‘Then for each local region, p,+ pau, 
remains constant on a C,. ‘This leads naturally to the suggestion that 
p,+ pau, varies very little along a C, even when large regions and finite 
changes are involved. If this were so we could argue that all C., character- 
istics start in x < 0 and p,+pau, = 0 there; hence, p,+pau, = 0 along each 
C , and in particular where they meet the shock. Roughly speaking, p,+ pau, 
would play the role of the Riemann invariant in simple wave theory. Indeed 
in a negative simple wave p,+ pau, is constant. ‘Thus, for example, in the 
simple wave formed by the C_ in x < 0 in figure 1, p,+ pau, vanishes, even 
though the total changes in p and u may not be smail. ‘This provides useful 
support for the contention that p,+pau, is small. However, in the general 
case the possibility that small changes may accumulate over a large region 
cannot be overlooked. Presumably, the result would only hold for flows 
which are slowly varying in some appropriate sense. Buta precise definition 
of ‘slowly-varying’ which would include cylindrical and spherical shocks, 
whose strengths ultimately become infinite, is far from clear. ‘The author 
has so far been unable to make further headway in this direction. ‘The 
expression for the rate of change of p,+ pau, on a C, is quadratic in the first 
derivatives of the flow quantities and so is of smaller order when these 
derivatives are small, but the expression is quite complicated and does not 
seem to suggest anything else significant. It is true that it could be used 
to find a correction term in going to a next higher approximation. But 
the smallness of the correction will only confirm the accuracy of the lowest 
order approximation without throwing further light on the reasons 
for it. 

This can be said of all the other approaches tried so far. However, 
one of them is quite systematic and gives a series of successive approximations 
with (9) as the first approximation. Asa mathematical procedure for finding 
the answer it is satisfactory, but it does not explain why the first approxima- 
tion isso good. ‘This is because it works essentially on the first factor in (12) 
instead of the much smaller second factor. ‘Thus one feels the method 
cannot be the right one. It is, however, closely related to Butler’s 
calculations for the converging cylindrical shocks and so is worth describing. 
‘This is done in the next section for those particular problems. 


3. CONVERGING CYLINDRICAL AND SPHERICAL SHOCKS 
We must first consider Guderley’s solution briefly. Guderley (1942) 
studied the limiting case of a very strong shock and found that there exists 
a similarity solution representing a converging cylindrical or spherical 
shock. For the cylindrical case we take equations (5) with A « (x) —.x) 
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corresponding to a wedge-shaped channel. Guderley’s similarity solution 
takes the form 
Vy X : (Xp ~ x r 


u ~~ Ff p = Py ~ —_ 


gE), p = poA(E), (17) 
where € = (x,—.x)/t* and « is an exponent to be determined. ‘The shock 
is a line € = constant = &,, say, so that its motion is given by the law 
Xy—x = &,t%. For example the shock velocity varies like (xy—.+)~" where 
n= (1—a)/x. When the expressions (17) are substituted in (5), ordinary 
differential equations are obtained for f, g, A and the solution satisfying the 
shock conditions at € = € is required. ‘lhe exponent « is determined in 
a rather strange way; there is only one value of « which gives a solution 
free of unrealistic singularities in the flow. For other values of «a fold occurs 














Reflected 
t shock 
x —- X 
fe) 

Limiting 
choracteristic 

/ incident 

shock 


Figure 2. ‘The (x, ¢)-plane for a converging cylindrical shock wave. 


in the (x,¢)-plane and the solution is no longer single-valued. This fold 
would be formed along a particular line € = €, which turns out to be the 
characteristic passing through the centre x»—x = 0 (see figure 2). To 
avoid the singularity the solution must have certain special properties at 
£= &, and these essentially determine x One of Butler’s suggestions 


$ >1 


(Butler 1954) for calculating the solution is to start with a series expansion 
near € = €,, thus ensuring regularity there, then the shock conditions at 
é €) are sufficient to complete the solution and determine «. 

We now consider this sort of approach directly in the physical plane 


without using the similarity solution. ‘This is independent of the assumption 
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of a strong shock and can be extended directly to other problems. Let 
p(x), u(x), etc., denote the values of flow quantities on the limiting 
characteristic. ‘Chen the characteristic is given by 7 = 0 where 


rat-f oS (18) 


_ yO + @®* 
Next, all the flow quantities are expanded in power series in 7: 
p(x, t) = pO(x) + tH (x) + 7°p(x) +..., (19) 
u(x,t) = u(x) + ru? (x) + r2u (x) +... ete. 


If we substitute these in the equations of motion and equate coethicients 
of 7 to zero, we get the characteristic relation 
dp du pa 0)244(0) l 


— +p'a’ — + =U 20 
' dx u™ +a (x9 —x) (") 





connecting the lowest order terms, and a set of equations for determining 
the higher order coefficients p'(x), p(x), etc., in terms of p™, p®, u. 
(This is in accordance with the general theory of characteristics.) Now, 
the boundary conditions (8) at the shock require 


) 





y—1 

(x) + T'(x)p™(x) +... = py a? {—— M2 - 7 —_\ 
p(x) + T(x)p(x) po) I 

pap y+1)M?* | 
(Of) (Lyf) a = ees? es! \ ? 
p(x) + T(x)p%(x) + = pu a5 + = (21) 
0( och, 7( :) (1)( a — - VM — | | 
u(x) 4 x)u™|(%)+... = a Tt a) | 


where 7(x) is the value of 7 at the shock, Le. 


x 1 1 ) 
io _ af dx. (22) 


For a first approximation, only the zero order terms, i.e. the values on the 
characteristic + = are retained on the left of (21). Clearly, then, on 
substitution into (20), we have exactly the simple rule of this paper. From 
this first approximation for p®, p, u, the coefficients p™, p', uw? can be 
found and so for a second approximation two terms can be retained on the 
left of each of the equations in (21). ‘This gives improved values of p™, 
p, uv in terms of M, which (on substitution in (20)) leads to an improved 
expression for M(x) and so on. Of course the improved function M(x) 
ditfers very little from the first approximation. It is an unexpected bonus 
in this approach. For, the whole approach is based on an expansion in 
powers of (w+a—U)/U (see (22)) which is not so very small. In fact 
the values of p™, p®, u change by the expected much larger amounts 
between the first and second approximations. For a strong shock, these 
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quantities are proportional to powers of (xy—.x), and for y = 1-4 the 
coefficients in the two approximations are : 





First Second 
approximation approximation 


| ane ‘ewer 

| p™ | ()-833 | 0-967 
ALD 6-000 7-420 

Th 0-833 0-784 | 





Putting it slightly differently, the first approximation takes just the values 
at the shock and the second approximation gives improved values at the 
characteristic. ‘The reason for the small change in the law of propagation 
for the shock is that the correction to the exponent is proportional to 


pP + pau 


and this is very small. In fact this quantity is the value of p,+ pau, on the 
characteristic! ‘Thus, the previous comments are vindicated. 


4. CONNECTION WITH SHOCK EXPANSION THEORY 

Shock expansion theory treats quite a different aspect of shock propagation 
but there are one or two points worth mentioning. ‘This method concerns 
propagation down a uniform tube when the effective piston motion varies. 
(Its practical importance lies in the applications to the analogous problem 
of supersonic aerofoils.) When the shock is weak, the flow is approximately 
a simple wave on the C,, characteristics with the Riemann invariant « given 
by 








and the entropy S equal to the undisturbed value S,. If the shock is not 
weak, a surprisingly good approximation to the pressures on the piston is 
obtained by assuming that the flow is a simple wave with 


2a 2a, 


ery aa es 





— Uy, S = §, 


where the subscripts indicate values just behind the shock initially (these 
are easily found from the shock relations and the initial velocity of the 
piston), ‘This is shock-expansion theory. ‘The formula actually used 
for the pressures on the piston is 

dp —pa du = 0, (23) 


because this holds for the simple wave. Since S = S, is exactly true on the 


piston, p and a can be expressed in terms of p, and (23) then determines 
the pressure variations in terms of the piston motion. Now, (23) is related 
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to the fact that the characteristic equation appropriate to the C_ is 
(Pr + (u—a@)p,; — patu,+ (u—a)u,} = 0, 
(for a uniform tube). ‘Thus if 
Pp; — pau, = 0, (24) 

then (23) follows. ‘This is the same condition as the one used in the work 
on non-uniform tubes but with the opposite sign. ‘The sign change is 
because the main disturbance is on the C,, characteristics here, but for the 
non-uniform problems (with a uniform piston motion) it is on the (¢ 
‘The same thing can be seen from Chester’s solution (14), (15), (16); we 
have 

Pr- PQ, = — (Uy — a, )G" (x — tu, — a3), (25) 
so that the reflected disturbance given by G is ignored in shock-expansion 
theory. In our discussion of non-uniform regions G is one of the main 
terms. We may also remark that shock- -expansion theory makes this 
additional approximation of neglecting G even in the small perturbation 
theory. The treatment of non-uniform regions is exact in that theory. 
For further details of shock-expansion theory see Eggars, Savin & Syvertson 
(1955) and Mahony (1955). 


KINEMATIC SHOCKS 


Analogous and much simpler shock problems arise in the subject called 
‘kinematic waves’ by Lighthill & Whitham (1955), and one might hope 
that discussion of these easier cases would give further insight into the 
more complicated cases. In fact it does not: the rule is correct but trivial. 
We consider, for example, kinematic waves in traffic flow on crowded roads. 
The density of the traffic, k, and the rate of flow across any section, g, must 
satisfy the continuity equation 


whe -. be = 0. (26) 


In addition qg will be some function of k, since the average speed of the cars 
(q/k) will fall as k increases. But if the road conditions are non-uniform, 
this relation will also include x. A shock wave, i.e. a discontinuous jump 
from conditions (9,,k,) to (qo, R2), travels with speed 
a ees - =r 

this is the only shock condition. ‘The problem is to find out how such a 
shock, initially travelling with constant speed along a uniform stretch of 
road, varies when it reaches a non-uniform section specified by 

= f(g»). (28) 
The solution is obvious: the shock continues to separate regions in which 
q takes the constant values q, and q, respectively. The corresponding values 
of k are functions of x given by (28), and (26) is satisfied. The shock velocity 
varies according to (27). 
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Now this solution is exactly the one given by the rule of this paper. 
‘The characteristic relation from (26) and (27) is 
f(g, x) ec ar 
: 4 ) 4 + vs i Or dq (), 
( q ( t COX 
Hence the rule takes g = constant along the shock, i.e. gy is constant. “Then 
since q, 1s given to be constant ahead of the shock, we have exactly the above 
solution. ‘The rule is correct for this type of wave motion but is so trivial 
that it does not help with the more complicated cases. 
FURTHER APPLICATIONS 
6. PROPAGATION THROUGH A STRATIFIED LAYER 
Suppose we have a stratified layer with the equilibrium values of pressure, 
density, etc. depending on vy. In general we can suppose that the layer is 
maintained by a body force F (which could also vary with x). Then the 
equations of motion are 


Pte up, + pu, si U0, 
1 m 
ut+uu,+-p,=F, } (29) 
p 


p,+up,—a@*(p,+up,) = 0. 
In equilibrium p = po(%), p = po(x) where 


he F. (30) 
Po dx 
The determination of p, in terms of pp is completed by giving the equilibrium 
entropy distribution. 

Chisnell (1955) considered this problem in the special case F = 0, 
Pp = constant (py not constant) by replacing the continuous density 
distribution by a piecewise constant function. The successive layers, 
with p constant in each, can remain in equilibrium because the pressure 
is the same ineach. But for varying fp it is not quite obvious how to proceed 
by that method, since jumps in fp, cannot really be allowed at the interfaces ; 
there is no force to balance the pressure difference (F is a body force). 
For the present rule, there is no difficulty. The appropriate characteristic 


relation is 


or 
a . 
dp + pa du - a dx = 0. (31) 
u+a 
The final step is to substitute the shock relations to determine the variation 
of \/J with x. In general the resulting first-order equation W(x) will require 


numerical integration. However, in the special case of very strong 
shocks, (8) simplifies to 
2 ' y+1 : een 
Pp —— py U®, poo —— py, u = — U, (32) 
y+l y—l y+1 
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> 


and we see that the third term in (31) can be neglected (since U is very 
large). ‘Thus the contribution of F is made indirectly through its control 
of the equilibrium density distribution py(v). With these approximations 
the law of propagation of the shock is 

U x py, pxp,”, 
where 





For y = 1-4, B = 0-2158. 
A specific example where this law may be useful is in finding the effect 

of the non-uniformities in the atmosphere on explosion blasts. In that 
case, F is the acceleration due to gravity, which is —g if x is measured 
upwards. If we consider an atmospheric layer in adiabatic equilibrium, 
i.e. Py = Kp, we have (from (30)): 
yKpy 


= €C — 2X, (34) 


where C is a constant. 


7. CONVERGING CYLINDRICAL SHOCK WAVES IN 
MAGNETOHYDRODYNAMICS 

This is a problem of some current interest which we can consider as 
a useful example where both changes in area and equilibrium density 
distribution occur. In such cases similarity solutions of Guderley’s type 
are limited to power law distributions of density (as well as to very strong 
shocks) but the methods of this paper can be applied to the general solution. 
The investigation will be based on the Lundquist equations of magneto- 
hydrodynamics (Lundquist 1952), and in particular the shock relations, 
etc., are all taken from a report by Friedrichs (1955). 

In terms of the velocity vector u, the magnetic field vector H, pressure p 
and density p, the equations are 


oH +Vx(Hxu) =0, (35) 

® +.¥. (pu) = 0, (36) 

— +(u.V)u Lup Hx (VxH) =0, (37) 
p p 

L +u.vp-a(2 -u. Vp) =0, (38) 


where yp is the permeability. ‘The Maxwell equation V.H = 0 is essentially 
included in (35) since the divergence of that equation gives 0(V.H)/ot = 0. 

For flows with cylindrical symmetry it is assumed that u is radial and 
that all flow quantities are functions of the radial distance r and the time f. 
The transverse and axial components H,, H. of H need not vanish, however. 
Indeed they are of primary importance since it is easily shown that the 
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radial component H, must vanish or the solution becomes trivial. For, 
if H, £ 0, the condition that the # and z components of H x (V x H) must 


vanish in (37) requires that rH, and H. be independent of r. In turn these 


lead to artificial forms for u and the only feasible case turns out to be 
H, = H,=0. But then the flow problem is independent of the magnetic 
field. ‘Thus, only the case H, = 0 is considered. 


With these conditions of symmetry, equations (35)—(38) reduce to 





OH, 0 
—§ + —(H,u) = 0, 
ot or 
ow. i@ 
— +-—(rH,u) = 0, 
ot ror ° 
0 Ou u 
cP u—t Pras = 0 
ct o or r 
1 Ou ] C ] L C - ms . H 
a iba So +o 24 H2)4° 0, 
ot cr po 2 ocr TS 


In the equilibrium state ahead of the shock, these equations reduce to 
- : ‘H* 5 
Pot Sp(H?, + H*) +p |\—dr = constant, (39) 
which determines the pressure distribution in terms of the magnetic field. 
(Subscripts zero are used for equilibrium values.) ‘The density distribution 
is deduced from py and the initial entropy distribution. In the usual case, 
the entropy will be uniform so that py & pz. 

‘Turning now to the shock relations, the general form simplifies in this 
case since the magnetic field has no normal component. ‘The appropriate 
relations giving the values behind the shock in terms of the shock velocity U/ 
and the equilibrium values are 

H(U—-—u) =H, U, 
p(U—u) = p, U, 
p+ SpH? + p(U —u)? = po + bu? + py U?, 
1([ 9 yp Ht? _ arte YP o wH* 
1 u)*? + —=— +4 -= $U24+ ——— + —’. 
(y l )p p (y —] )Po Po 
It is convenient to be able to solve these relations and express all the 
quantities in terms of a single parameter characterizing the shock strength. 
This can be done in terms of the parameter € = p/p). Then, 








¢ 
‘ be ae oe 
p Pos a = H,é, w= Us ot 


U2 a | 24 2 (1 Y\ ¢ | (40) 
_ ns | (TH -- [)* < ~-—~J&+ = ff, : 
(yv+1)-(y-1DEL" * | 2 





2p)(€—1) [ Me. ak oe | 
pee... at) Ae ek! 
P=Po* C+ 1)-@-IeL” 4 | 
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where ay is the sound speed \/(yfo/p)) and dy is the Alfvén speed 1/(4H?/py) 
(i.e. the speed of magnetic waves in an incompressible fluid). ‘There are 
two possibilities leading to strong shocks, 1.e. large values of p/p); either 
€ is close to the value (y+1)/(y—1) or 6?/a? is large. The former is the 
usual case found in gas dynamics, but it is an interesting fact that in 
magnetohydrodynamics strong shocks may also arise when the magnetic 
field is very strong for any value of €> 1. When & + (y+1)/(y—1), the 
shock relations take the simple form 








— ed 

p=: 7 H =* H,, 

) y- : y-1 41 
20 172 2 ou . sai 

p~ y+l ’ u~ y+i 


It should be noted that in this limiting case, the expressions for p, p, u, 
in terms of U are independent of the magnetic field. 

Now we apply the rule described in the earlier sections to determine the 
motion of a converging shock. ‘The characteristic equations corresponding 
to the first pair in (6) are 


10 @ 12a. a... 
(ai ; =) (0 smHG + 3uld:) + oe(—— st =)e+ 


_peul pHyuye 











O (42) 
ut+er r ute 
where 
sappme ll ee. (43) 
pp 


In this section the particle velocity u is taken as positive in the direction 
of increasing r.. For a converging shock, therefore, u is negative and we 
must use the relation for a negative characteristic to determine the shock, 
1.2; 


dp +uH,dH, +uH.dH.—pe du { :! fe faa - a - =(0. (44) 


The final step is to substitute the shock conditions (40) into this relation. 
A first-order equation for €(r) is obtained, and this determines the shock. 

In general, the equation will require numerical integration. But in the 
special case of very strong shocks with € > (y+1),(y—1), the simplified 
shock conditions noted in (41) may be used. From these conditions we 
see that the corresponding approximations in (44) reduce it to 


pau dr 


(a Tr 


dp—pa du+ F = Q, 





since the other terms involving the magnetic field are of smaller order in 
the shock velocity LU. Now, when (41) is substituted in this equation, the 
resulting equation for U(r) integrates immediately to 

U x px *r-r, (45) 


F.M. Z 
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where 
a 7 es ‘, 
Bsa ds |(4) - ¢ eeauerts (2) . £24 |(-4) 
. ez, L Ae i ak 2 
(46) 
‘The corresponding law for the pressure behind the shock is 
pa p! si) Sia (+7) 


‘The exponent # is exactly the one appearing in the solution for a converging 
cylindrical shock in a uniform medium (see equation (11)) and § is the 
exponent that was found in (33) for a plane shock in a non-uniform medium ; 
(45) shows the combined effect. It 1s interesting to note that for a strong 
shock the magnetic field enters the law of propagation only indirectly through 
its control of the equilibrium density distribution p,. This is exactly 
analogous to the situation found in the last section for the effect of the body 
force F. No doubt (45) and (47) apply in general to any case when the 
distribution is not uniform whatever mechanism maintains the density 
distribution. But for weaker shocks, the details of the particular force field 
would be required. Since 1—28 > 0, a general consequence of (47) is that 
the shock will be relatively strengthened or weakened according as the 
density p, increases or decreases towards the centre. 

The law of propagation of the shock wave can be given explicitly in the 
case of weak shock waves also. ‘This case is less interesting and only the 


final result is noted. It is found that 


P—-Po " 1/2.~3/2,— 12 
Po 


&. SHEAR LAYERS IN SUPERSONIC FLOW 


‘The problems in steady supersonic flow are analogous in many ways 
to the unsteady problems considered so far. ‘There is, however, one important 
question that arises in steady flow problems; the flow downstream of the 
shock may be either supersonic or subsonic. “The questions and methods 
discussed in this paper apply where the flow is supersonic; if the flow 
downstream of the shock is subsonic the methods are meaningless. With 
this restriction the problems can be solved in the same way and the 
perturbation methods etc. discussed earlier go through and give similar 
results. Accordingly the discussion will be brief. 

In this section we note how Moeckel’s original results can be obtained 
by application of the characteristic relation. We consider the plane flow 
represented in figure 3. The incident flow is uniform in y < 0 but varies 
with y in y > 0. The only requirement is that the pressure py be constant ; 
then, the density distribution p)(v) and velocity distribution gq (v) can be 
arbitrary. (In reality, of course, viscous forces will arise in such a shear 


layer but we must confine attention to the case when these can be neglected. ) 
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In supersonic flow the derivation of the characteristic relations from the 
usual equations of motion involves much more manipulation than in the 
unsteady flow discussed in §3. The details are given in most books on 
compressible flow (aithough not all books give the general non-isentropic 











Figure 3. Interaction of an oblique shock with a shear layer 


form required here) and it will be sufficient to quote the result. On a 
positive characteristic (C, in figure 3), the flow quantities must satisfy 
the relation (see, for example, Howarth (1953)) 
ld dé dy 
ithe 0 + ——_——_ = J). on = = tan(p+9), (48) 
yp  cospsing dx 
where p is the pressure, 4 is the angle of the velocity vector to the x-axis, 
wis the Mach angle sin~!(a/q) where a is the sound speed \/(yp/p). ‘The 
rule is to apply this relation between p and @ along the shock. 
The shock conditions may be written 











2v, 4 oe 1 
p = Pos ——~sin? B— ¢ Fe 
ly+ ide) 
y-1 2 a. ti 
ee ( i ee TT , (+9) 
y-1 e ps 1 | 
¢ f — rn 2 I Nar Ey 
tan(f —6@) oe ane, y+1M?sinf cos f 


gcos(B—@) = qycos f, 


where § is the angle of inclination of the shock to the x-axis. From (49) 
all the flow quantities can be expressed in terms of § and the upstream flow. 
The shape of the shock is found by substituting these expressions into (48). 
A first-order equation is obtained for 8 as a function of y. Of course, any 
Z2 
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of the flow quantities can be chosen as the variable instead of 8 and a similar 
equation obtained. 
Moeckel introduces the following notation for the solutions of the 
relations (49) 
P= Pofi(“ B), 
@=flM),B), | (50) 
y(sin «cos 1) — fs(Mp, B). 


When these are substituted in (48) (noting that py is constant), we have 


dp l chy f Che ( I ch ' Cfe\ 4 
a - — — +f,.— -— +Ja—]. 1 
dM, (; 0M, Is: =) fi &p Is F 3) (51) 


I'his is Moeckel’s result. 


9, AXISYMMETRICAL DUCT IN SUPERSONIC FLOW 
\n interesting application in supersonic flow is to the shock formed 
in the entry of an axisymmetrical duct. If r measures distance from the 


ixis, the appropriate characteristic relation is now 


—— — = 0. (52) 


yp COS 4 SIN pL cos psin(u—§) r 





‘The change in sign from (48) arises because the shock instde the duct is 
formed by the negative characteristics; since @ is negative in this case the 
signs of (48) and (52) eventually agree and the only essential difference is the 
additional term in (52) due to the cylindrical geometry. The shock shape 
is then obtained by substituting the shock conditions (49) into (52). In 
this case we assume uniform flow upstream of the shock so that M, etc. 
are constants. ‘Then, the equation for f(r) takes the form 
dp 


r = = —f(), (53) 


where f(8) is a known function. It is found that f(8) > 0 so that £8 increases 
as r decreases towards the axis; hence the strength of the shock increases 
towards the axis. However at some point before the axis is reached the 
flow behind the shock becomes subsonic and certainly (53) does not apply 
beyond this point; in fact, f(8) becomes imaginary. In the real situation, 
Mach reflection occurs and the flow pattern becomes very complicated. 
It is easy to see that the Mach reflection must occur somewhere before 
this ‘sonic point’ is reached because the reflected shock must have 
supersonic flow upstream of it. However it seems reasonable to assume 
that (53) applies until the triple point is reached. 
If we denote the radius of the duct by r, and the initial angle of the 
shock by fp, (53) gives 
log 2 =| f(8) dB. (54) 


ar z 


The initial shock angle 8, is a function of the initial angle 6, of the wall 
of the duct (it is given by the third equation in (49)). Thus, according 
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to (54), only the initial angle of the wall duct affects the shape of the shock 
appreciably. This corresponds to the result for cylindrical shocks that the 
motion of the shock is nearly independent of the piston motion (see § 1). 
At the sonic point, 8 takes a value 8, depending only on M, and y, and the 
position of this point is given by r = r, where 


r “Bs mn 
log =| f(8) dp. (5 
‘, 


~ Bo( Ao) 


wn 


) 


Asa check on the present method, the predictions of (55) can be compared 
with calculations of the shock shape by the numerical method of character- 
istics (Ferri 1946). In particular, it is convenient to check the dependence 
of r)/r, on 4. The integral in (55) was calculated for a few values of 4, 
using VW, = 1-6, the value chosen in Ferri’s calculations, and the results 


are given in the following table. 





Gy 14-24 11-62 


Filho 1-000 0-659 | 0-337 0-048 








A graph of the function is given in Ferri’s paper and the values given above 
fit almost on the curve; the curve is shown in figure 4+. This problem is 
a fairly severe test of the simple rule presented here and the very close 
agreement with accurate calculations gives ample support of the practical 
usefulness of this method. 
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Figure 4. The variation of flow angle @, with distance 7, from the axis in supersonic 
flow inside a duct. 


10. PROPAGATION OF BORES IN SHALLOW WATER 
The application to this case is quite straightforward but the results 
require special comment. 
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If hA(x,t) denotes the depth of the water, Ao(x) the undisturbed value 
ahead of the bore, and u(x,t) the particle velocity, the equations of the 
shallow-water theory are 

mt+i(hy+n)u}, = 9, 


(36) 
u,+uu,+gn, = 0, 
where » =h—h,. The appropriate characteristic relation is 
go dh, 1 
du t 2dc Rat ’ (D/ ) 
7-3 
where ¢ vy (gh). ‘Vhe bore conditions are 
oh(h, +h) a 
[ = (f (My +") ‘ ( SS ) 
2h, 
h—-h, ;,; se 
an ; (59) 
h 


It is convenient to introduce VW = U/, (gh) and determine first the variation 
of with x. Then the height of the bore is given by 

n = h—h, = 2h,(M?—-1) (60) 
ind the bore velocity by 











: (gh) = M /(2M?-1). (61) 
From (59), (60), and (61), 
Tah} = ,/(2M?-1), 
(62) 
u 2M(M?-1) 
v(ghh)  yv(2M?-1) 
Substituting these expressions in (57), we obtain 
1 dhy 2(2M* + 2M?—2M—1)(4M*+ 4M3—3M?-2M +1) (63) 
h,dM (M?—1)(2M*—1)(2M*+6M?+2M?-—3M-—2) ; 


The range of variation of M is 1 < M < » and in this range all the 
factors in (63) are positive. Therefore, WM always increases as h, decreases 
and vice versa. For weak bores (n/4y < 1), M is near 1 so the asymptotic 
form for (63) is 


hydM 


iy * 3 
5M-1- 

Hence, 

M-1« h,*4, noc a, (64) 
This agrees with the result of the linear theory. For strong bores 
(y hy > 1), M is large and the approximate form of (63) is 

1 dhy _ _* 

h@dM@ M 


Hence 


Moa h~™4, nx hi?, (65) 
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It should be noted that for weak bores 7 increases as fy decreases, but 
for strong bores » decreases as fy decreases. ‘The critical value MW, of M 
which separates these two regimes is the value for which 

dy dM 


1 _ 9M2—1)+4h. M&™ = 0. 
5 eee 


‘The equation for .W/, is then found from (63) and the value for M, can be 
calculated. It is found that V7, is approximately 1-2; the corresponding 
value of 7/Ay is 0-9. 

This result for the variation of » with Ay has strange consequences 
when we apply the above theory to consideration of a bore coming into the 
shore line of a sloping beach. Assuming that the bore is sufficiently weak 
initially, it will increase in height as expected, but eventually, although 
its strength 7/h) continues to increase, its height » will decrease. In the 
final stage, » x h'* as given by (65) so the height of the bore tends to zero. 
This is certainly not in accord with preconceived notions of what should 
happen. We might consider that the lowest order approximation is not ade- 
quate to treat the extreme case in which Ay -> 0. It would not be surprising; 
in the analogous problem of a shock moving into a stratified medium, Chisnell 
found that it was necessary to go to a second approximation (i.e. include 
‘re-reflected waves’ in his theory) in such extreme cases. However, for 
this water-wave problem the idealizations of the basic shallow-water theory 
and bore conditions are drastic ones, and we can see immediately that they 
lead to trouble. For, one expects that in the true situation the bore 
velocity U and height h approach finite values at the shore. But, if this 
is so (58) cannot hold as A, -- 0. Thus any theory which includes (58) as 
a basic condition cannot lead to the expected results. In fact (65) does 
predict a finite value for U so (63) gets one of the expected results. For 7 
to have a finite value as h, — 0, U must become infinite like 4, !? according 
to (58); this result would have been unsatisfactory also. It is easy to see 
why shallow-water theory may break down when strong bores are involved. 
In shallow-water theory vertical accelerations are neglected and the motion 
is assumed to be essentially horizontal. It is possible, therefore, that the 
theory breaks down when the large change in water levels at a very strong 
bore is considered. 

On the other hand, even though the idealized model may not describe 
reality very well it does conserve mass, and the feeling that the height of the 
bore should be finite is based partly on this. For the height of the water 
surface must be increasing on the whole, since we are assuming that there 
is a region of uniform depth away from the shore providing a continual 
inflow of water. If the bore height tends to zero it must be followed by 
large (continuous) increases 1n depth. 

To settle these points much more extensive investigations are necessary. 

This paper represents results obtained at the Institute of Mathematical 
Sciences, New York University, under the sponsorship of Contracts 
N6ori-201, T.0.1 and Nonr-285(06) with the Office of Naval Research. 
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The effects of radiative transfer on turbulent flow 
of a stratified fluid 


By A. A. TOWNSEND 
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SUMMARY 


Assuming local thermodynamic equilibrium in the fluid, an 
expression is derived for the rate of destruction of the mean 
square of the temperature fluctuations by radiative transfer of 
heat. This takes a particularly simple form (a) if the fluid is 
effectively transparent over distances equal to the scale of the 
turbulent motion, when the effect appears as a decay time for 
temperature fluctuations from the mean, and (d) if the fluid is 
effectively opaque, when the effect is of an increased conductivity 
due to radiation. A theory of the interaction of the temperature 
and velocity fields developed in a previous paper shows that, if 
the radiative effects are relatively weak, a sudden collapse of the 
turbulent motion occurs while the flux Richardson number is 
still less than one. If the radiative effects are strong, the turbulent 
intensity approaches zero as the flux Richardson number 
approaches one. ‘The effects of radiation are always to increase 
the critical value of the ordinary Richardson number. Criteria 
for fully turbulent motion of an unrestricted flow are given in 
terms of the gradients of mean velocity and mean temperature 
and of the rate of radiative cooling. ‘The relevance of these 
calculations to motions of the atmosphere is briefly discussed. 


1. INTRODUCTION 

In the past decade, evidence that the upper atmosphere between seventy 
and one hundred kilometres above sea-level may be in turbulent motion 
has been accumulating. Most of the evidence is indirect and based on the 
irregular fading of radio echoes from the ionosphere (e.g. Briggs, Phillips 
& Shinn 1950), but Liller & Whipple (1954, Vol. 1, p. 112) have made 
observations of luminous meteor trails which show the existence of irregular 
velocity gradients resembling closely those found in turbulent flows. The 
nearly constant composition of the atmosphere from sea-level to these 
heights also shows that mixing on the molar scale must occur sufficiently 
frequently to counteract diffusive separation of the atmospheric components, 
and these mixing motions might be turbulent motions. 
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As most of this evidence of turbulent motion is indirect, it is useful to 
ask whether our present knowledge of turbulence would lead us to expect 
turbulent flow with the gradients of mean velocity and temperature occurring 
in this region. ‘The criterion for turbulent motion of air of sea-level density 
is usually expressed as a critical value of the Richardson number and, 
although there is considerable disagreement over the critical value, a simple 
ipplication of this criterign to the upper atmosphere makes turbulent 
motion appear very unlikely*. However, the derivation of this criterion 
neglects radiative transfer of heat, a process which is known to dominate 
the heat balance in the upper atmosphere, and its inclusion will cause a 
reduction in the magnitude of the buoyancy forces which are responsible 
for the inhibition of turbulent motion. In general, radiative transfer tends 
to destroy the effects of density stratification caused by temperature gradients 
as Goody (1956) has shown in his investigation of cellular convection 
between parallel horizontal planes. ‘The purpose of this paper is to discuss 
the turbulent motion in a stratified fluid with appreciable radiative transfer 
of heat using a similar analysis to that used by the author in a previous 
discussion of stratified flow with negligible radiative transfer, and to derive 
criteria for the maintenance of turbulent motion. 


2. ‘(URBULENT TRANSFER RATES IN STRATIFIED FLOW 

In a recent paper (‘Townsend 1958), the interactions between the fields 
of velocity and temperature in a stratified flow were considered, and relations 
were found between the levels of the turbulent fluctuations and transport 
rates and the gradients of mean velocity and mean temperature. ‘These 
relations form the basis of the present paper but, before quoting them, it 
may be useful to set out briefly the steps in their derivation. 

We consider the steady flow of a nearly perfect gas, for example air 
containing only small proportions of water-vapour, carbon dioxide, ozone, 
etc., with velocity variations small compared with the local velocity of sound, 
with a length scale small compared with the scale height of the atmosphere, 
with temperature variations small compared with the absolute temperature, 
and unaffected by the rotation of the earth. Without serious loss of 
generality, we may suppose the mean flow to be nearly unidirectional and 
horizontal with an appreciable gradient only in the vertical direction, e.g. a 
horizontal mixing layer between two streams, and describe the flow in the 
usual coordinate system with Ox in the horizontal direction of flow and 
with Oz vertically upward. Then, 

U', Ware the components of the mean velocity parallel to Ox, Oz 

respectively, 

u,v, w are the components of the velocity fluctuation parallel to Ox, 

Oy, Oz respectively, 
T is the mean absolute temperature, 


* For example, for a zero lapse rate and a wind shear of 10 metre sec-! km7@! 


the Richardson number is 3:5. 
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# is the local temperature fluctuation, 
v is the kinematic viscosity, 
« is the thermometric conductivity (heat diffusivity), 
pis the mean density, 
c, 1s the specific heat at constant pressure, 
g is the downward acceleration due to gravity. 

To the approximation implied by the restrictions set out above, the 
equation for the kinetic energy of the turbulent velocity fluctuations* is 
aig) aU , ake) , AGP) , 2 (Ia, 1\_ 2 

p—— + U —=— + W -— G+ = (5 Fe: ~ Pi) = 2 dw-—e, 

ct Cz OX os ( 2 p 1 
(2.1) 
where g? = u® +77 +w? and «€ is the rate at which turbulent kinetic energy is 
being converted to heat by the action of viscosity. A second equation of 
physical importance is the equation for the mean square of the temperature 








fluctuation. It is 


O( Ly — /oT g _0( 16" , 0( 462 ( a ——— —s 
eh, ) - Aw | ae a ca al ) + sa ll + — ($6?w) = K6V70— Be, 
ct ” ON C3 os 

(2.2) 
where — £6? = #6(pc,,)~! and # is the net rate of gain of heat per unit 
volume by radiative exchange. ‘These two equations are, respectively, 
direct deductions from the equations of motion and from the equation 


for the internal energy (heat equation). 

To simplify the form of these equations, we make use of the observation 
that, in free turbulent flows remote from solid boundaries, there is a sharply 
defined surface separating turbulent fluid from the surrounding undisturbed 
fluid and that the turbulent fluid is remarkably homogeneous in intensity 
and scale over any section of the flow (‘Townsend 1956). ‘This means 
that there is a physical meaning in using averaged values of the turbulent 
intensity and of the mean square temperature fluctuation, defined as 


*- 20 -o 


I = l = 
— 2dz and — 6? dz, 
D}_.! D)_, 
where D is the mean vertical extent of the fully turbulent fluid. If 
equations (2.1) and (2.2) are averaged in this way, they become 











oe al ] 0 ~s - a 
uw — + — (4Uq*) = 2, dw-e (2.3) 
oz ox” : 7 
and 
= fcr . - 2 ps5 3 na 
6w{— +2) + — (4U6?) = KOV70—- Be, (2.4) 
ez & Ox 


where all the quantities are to be understood as suitably averaged values. 


* In a steady flow, the time derivatives of equations (2.1) and (2.2) are zero. They 
are included to make clearer the meaning of the equations. 
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We also neglect the terms representing gain or loss through advection, 
i.e. 0/0x(} Ug?) and ¢/éx(}U6?), and arrive at the comparatively simple 











equations, 
a ae 
uw — = 2 Ow-e, (2.5) 
Cs 7 
afer. 2 ACI __ Rat 
6w{— +2) = «OV70— Be". (2.6) 
OB ly 


It will be noticed that the radiative properties of the fluid do not enter 
directly into the equations for the turbulent energy and the convective heat 
transfer enters only through the buoyancy term (g T)dw. Itisa plausible 
assumption that the general nature of the turbulent motion, i.e. its 
geometrical properties but not its intensity or scale, depends only on the 
relative magnitude of the buoyancy term, that is, on the flux Richardson 


number 





a eng ol 7 
R, = ap Ow [ uw —, (Z./) 


which is the ratio of the rate of loss of energy through working against 
buoyancy forces to the generation by working against Reynolds stresses. 
Since energy dissipation is essentially positive, the flux Richardson number 
must always be less than one. 

‘The equation for the mean square temperature fluctuation (2.6) 1s a 
relation between the velocity field and the temperature field produced by 
its interaction with a gradient of mean temperature. Appreciable radiative 
heat transfer leads to added destruction of temperature fluctuations from the 
mean, to a lower mean square temperature fluctuation and to a lower absolute 
value of the convective heat transfer, #z. Whether the gradient of potential 
temperature is positive or negative (stable or unstable to convective 
disturbances), the effect of radiative heat transfer is to reduce convective 
heat transfer in a given velocity field and mean temperature gradient and 
so to reduce the contribution (negative or positive) to the kinetic energy 
of the motion. If the flow is convectively unstable, radiative transfer reduces 
the instability as Goody (1956) found in his investigation of the effect of 
radiative transfer on convective instability between parallel horizontal 
planes. If the flow is stable, radiative transfer reduces the stability. Indeed, 
if radiative transfer were infinitely rapid, no motion of the fluid could cause 
any departure of local temperature from the mean value appropriate to the 
position in the flow and no buoyancy forces could be generated. 

Equations (2.5) and (2.6) contain terms representing the rate at which 
viscosity converts turbulent kinetic energy to heat and the rate at which 
conductivity destroys temperature fluctuations from the mean. In fully 
turbulent flows, these rates are known to be independent of the actual values 
of the viscosity and the conductivity, and are determined by the large-scale 
characteristics of the turbulent motion, that is, by the intensity and scale 
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of the motion (Townsend 1956). This may be expressed by writing these 
rates as 

e = (w*)P2L-} (2.8) 
and —K0V70 = 167(w?)!?L751, (2.9) 
where L., L, are defined by these equations but are known to be nearly 
equal to the integral scale of the turbulence*. 

It is now possible to obtain solutions of equations (2.5) and (2.6), 
expressing the convective heat transport and the Reynolds stress in terms 
of the gradients of mean velocity and mean temperature, the logarithmic 
cooling rate by radiation f, the dissipation lengths L, and L,, and the two 
correlation factors, 


| ee (2.10) 





ed ae (2.11) 


The most useful form relates the ordinary Richardson number 





a (oT g\ //oU 
Dg S46 2 2.12 
i f i (5 £\, ( og ) ( ) 
to the flux Richardson number, and is 

H L, k2R,\12 
sas one = _19 6'*@ i 9) 4 
R,= 5 I (1 12 ea) | (2.13) 
wdhcke Pe (2.14) 


k, L. oU/dez 
It should be emphasized that the validity of equation (2.13) depends 
on only one assumption, that the effects of mean flow advection can be 
ignored, and that the quantities occurring in it are not local values but 
values averaged over a whole section of the flow. Its physical importance 
depends on the meaning given to these average values by the effective 
homogeneity of the flow and on the making of plausible assumptions about 
the variations of the non-dimensional ratios, k,,, k, and L./L,, with stability. 
3. CONDITIONS FOR THE MAINTENANCE OF TURBULENT MOTION 
From equation (2.13), we see that real values of the flux Richardson 
number are only possible if 
ea H?k,, L, (3.1) 
i= 2... «© J. 
12k, L, 
* The factor 4 in (2.9) appears since, in any flow, the ratios of the rates of destruction 
to the respective intensities are roughly equal, i.e. 
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Assuming that k,, k, and L/L, do not vary greatly with stability of the flow 
(reasons for supposing this in fud/y turbulent flow are given in ‘Townsend 
1958), this is equivalent to 


R, < 3H. (3.2) 


These are limits set by non-existence of physical solutions of equations (2.5) 
and (2.6) for Richardson numbers which do not satisfy them, and are 
additional to the limit set by the energy equation alone, 


R, <1. (3.3) 


For weak radiative transfer (H < 2, 8B < !k,(L./L,) cU/cz ), the limit 
expressed by (3.1) and (3.2) applies, while for strong radiative transfer, 
the limit is given by equation (3.3). 

The physical situation that leads to a double restriction on the possibility 
of turbulent flow may be made clearer by considering the equation for the 


turbulent kinetic energy 


—ocl g— (q2)32 - 
w= = 6+ =. (3.4) 
uw = ; L 


and, in particular, the variation of the terms with assumed turbulent 
intensity for given gradients of mean velocity and mean temperature. 
An intermediate step in the algebra leading to the relation between the 
two forms of the Richardson number is (Townsend 1953) 








F —|0U 3(g/T)k2 L,(w?)!?(eT/cz+g/c,) (w?)3!2 (3.6) 
Rw | = - +4 ; ae 
ut Oz 1 3 BL ,(w*) 1/2 L. 


Consider now the way in which the two terms on the right, representing 
work done against buoyancy forces and loss of energy by dissipative processes, 
vary with turbulent intensity. ‘This is done most conveniently by considering 
the ratio of their sum to the term on the left, the rate of production of 
turbulent energy by working against the Reynolds stresses. In figure 1, 
the variation of these ratios with (w?)!?/k, L,.0U/éz ) is shown for negligible 
radiative transfer and various values of the Richardson number, R;. If a 
stable value of the turbulent intensity exists, it must satisfy the energy 
equation, i.e. the sum of the ratios must equal one, and the total rate of 
energy loss must increase with turbulent intensity. It is clear from the 
diagram (i) that the sum of the ratios always has a minimum value, (ii) that 
the minimum value will exceed one in strongly stable flows, (iii) that the 
turbulent intensity does not become zero as the critical condition is 
approached, and (iv) that, in the critical flow, the flux Richardson number 
(which is the ratio of the work done against buoyancy forces to the energy 
production by shear) is less than one. ‘This behaviour is typical of flows 
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with weak radiative transfer and the condition (3.1) for turbulent flow is 
simply the condition that the minimum value of the sum of the ratios should 
be one or less. 

In figure 2, the sum of the ratios is shown as a function of turbulent 
intensity for critical conditions and increasing radiative transfer. The 
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Figure 1. Variation with assumed turbulent intensity of the ratio of the total rate 
of energy !oss to the rate of gain from the mean flow (no radiative transfer). 
(The numbers opposite the curves are values of (Lg k3/L, k2)R;, which in 
the text is equated with the Richardson number. The marked points indicate 
the stable configuration (if any) of the turbulent motion for the Richardson 
number concerned.) 


position of the minimum moves to smaller values of the intensity and 
becomes zero when § = 3k,(L,/L,) 0U/ez. For more intense transfer, 
the sum of the ratios always increases with turbulent intensity and the 
minimum possible value occurs at zero intensity. 
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The alternative conditions correspond with the suppression of the 
turbulent motion in two distinct ways. In conditions of weak radiative 
transfer, loss of energy through buoyancy forces is the dominant dissipative 
process for low values of the turbulent intensity as ordinary turbulent 
dissipation is for high values. In the critical condition, the intensity is 
neither so low that buoyancy forces could destroy the motion nor so high 


—_—____— 9 —_—___ 
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Figure 2. Variation with assumed turbulent intensity of the ratio of the total rate 
of energy loss to the rate of gain from the mean flow for the critical condition 
with various amounts of radiative transfer. (The marked points indicate the 
configuration and intensity of the turbulence in the critical condition for the 
various values of the radiation parameter H.) 


that turbulent dissipation would, but has an intermediate value. As the 
limit is passed a sudden collapse of the turbulent motion will occur. Since 
the effect of radiative transfer is to destroy temperature fluctuations and to 
reduce buoyancy effects, the loss of energy through buoyancy forces is 


reduced, particularly at low intensities, and for sufficiently intense radiative 
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transfer the critical condition has zero intensity. This corresponds with 
a condition of ‘just no turbulence’. 

The criterion for the maintenance of turbulent motion thus takes two 
forms, depending on the magnitude of H. The first is 


PRL, 


crit. ~ 12, 


(R,)., =H, (R) 


for H <2 or 


= 9) 
= 
~ 
K 





joUjoz|) 3 Le * 
The second, appropriate to conditions of strong radiative transfer, is 
k, £6 
(Ry erie = 1, (Ry a 


crit. ~ B2 9U oz 


\| 


for H > 2 or 


8 oe 


ae ot Se = ms 
dUjesi 3 LL, ° 





Obviously the Richardson number is not a convenient parameter for 
describing flows with strong radiative transfer, and these last conditions 
are better expressed as 





eam cs i) 
Pt ie 
és | crit. 


for 








4+. ‘THE EFFECT OF RADIATIVE HEAT TRANSFER ON 
THE TEMPERATURE FLUCTUATIONS 
In the previous sections, the quantity H (defined by equation (2.11)) 
appears prominently and it is a measure of the ratio of the logarithmic 


rate of cooling by radiation 8 (equation (2.4)) to the mean rate of shear. 


To compute 8, we need to know something about the dependence of -7, 
the net rate of heat gain per unit volume by radiative transfer, on the 
temperature distribution in the fluid. In general, # is a complicated 
expression depending on the difference between the temperature radiation 
from the fluid volume and its rate of absorption of radiation from other 
parts of the fluid and from outside the fluid, but our present concern is only 
with fluctuations of temperature and so we do not need to know the nature 
or quantity of the constant radiation from outside the flow. ‘To calculate 
the radiation from other parts of the fluid, we assume that the gas is 
everywhere close to thermodynamic equilibrium and that we may neglect 
the effects of unmodified scattering. 


F.M. 


Ny 
- 
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The equation of radiative transfer is (Chandrasekhar 1950) 

1 di, 
kp dr 
where k, is the mass absorption coefficient for radiation of frequency », 
I, is the specific intensity, i.e. the flux of radiant energy per unit solid angle 
per unit frequency interval, 7 is distance along the direction of propagation, 
and B,(7) is the Planck function for black-body emission. ‘This equation 
may be integrated to give the intensity at a point as 





£ 


I,= | B,(T,)e-*v’'k, p dr, (4.2) 


where r is the distance between the point concerned and the moving point 
of integration, and the density is assumed to be substantially constant. 
The total intensity from all parts of the fluid is then 


~~ 
+ 
o>) 
~— 


| I, dw = | B(T,Je kor ®y 5 dV (x), 
: i 


where the integrals extend respectively over all solid angles and all the 
fluid, the subscript 7 denotes values at the element of integration, and J’(r) 
is an element of volume at the position r. From this it follows that 
absorption of temperature radiation from the surrounding fluid proceeds 


at a rate per unit mass of 
, aa 1 : 
Bos B(T.)R,e K yer =p dl(r)dv, 
J di 


or, more concisely, of 


A Sette) 96 t ary), (4.4) 
7. e(pr)- fi 
where 
E(T,pr)= | BAT)l—e-*r) dv || BT) d 


is the emissivity of the fluid. Heat radiation from an element of temperature 
T proceeds at a rate 


4h(T)oT?, (4.5) 


KT) =| k, B(T)dv/ | BAT) dv Bow , 


C(pr) 


where 


Both E(7, pr) and k(T) depend on temperature, but if the larger values 
of Rk, are distributed fairly evenly within the range of the Planck function, 
B,(T), both will be slowly varying functions of .temperature. 

Neglecting the variations of E(7,pr) and k(7) with temperature, the 
net rate of temperature rise by internal radiative transfer is 


-_ are ( ((T rn o-— (T ' A)*] E (pr) 





p dI(r), (4.6) 


C A, re 


B(T), (4.1): 
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where dashes denote differentiation of E(pr) with respect to pr (note that 
E" (pr) is essentially negative). 

If the expression (4.6) is multiplied by the temperature fluctuation and 
the mean values taken, the result is the radiation term in equation (2.4) 
which represents the rate of destruction of 592 by radiative transfer. 
Assuming that @ < 7, this gives 





= OP as o E" (pr . ‘ 
= a ee Se, | (7, +6,)'-— (T+ 6)"]6 e () Fill (r). (4.7) 
pe, My . ‘= 
It is probable that the temperature covariance 6(x)#(x+r) is nearly 
symmetrical about a plane perpendicular to the gradient of mean temperature, 
and then 


eT | ae og 
— —— | (06,- 6] 








_ pm = =e") 5 aV(r), (4.8) 


7C 


p 
indicating that the destruction of }@? by radiative transfer takes place at 
a rate whose maximum value in a fluid of given properties is 

4o0T* — ( E"(pr) 
-——F | hal pdV(r) = 16 
p . - Cp 
This maximum value occurs when the scale of the temperature fluctuations - 


a 
k pe. (4.9) 





7C 


is so small that 96, becomes negligible before the transmission is sensibly 


different from one*. In general, the rate of destruction of 36? is 








saa —¢ i, = 
Be? = 16Fk a (4.10) 
a 
where : l ( 66.\ E"(pr) ; 
P= — — — — ) ——pdV(r 
47 } =) kr? , (*) 
1 (66, 16), E (er) 
1+ pdV 4.11 
Defining the average value of 66, ane over a spherical shell of radius r as 
06, 
or) = [ ase) (4.12) 
a 4nr* ! @ 


(dS(r) is an area element of the spherical shell), we find that 
F=1+ { Ov) ae pil (4.13) 


showing how F depends on the relative extents of the temperature corre- 
lation function O(r) and the transmission coefficient II (pr) = E"(pr)/E" (0). 
For two limiting conditions, the factor F may be expressed more simply. 
(a) If the fluid is effectively transparent, i.e. if 
-¢ 1 i: c x 
| Il(pr)dt=-—--=—> | O(r)d, 
pE"(0)~ Jo 
*In this event, the fluid is effectively transparent, radiation from other parts of 
the fluid is not absorbed, and equation (4.9) is a direct consequence of (4.5). 
2A2 
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then 


y) i w . 
pe) |" ovr) ar. (4.14) 


v 


(6) If the fluid is nearly opaque, 1.e. if 
| Il(pr)dr < | O(r)dr, 


only small values of r for which O(r) is nearly one contribute to the integral. 


cH oe Son 
O(r) =1- Ar? (=) 16, 


For these values, 





OX] | 
and so 
pa. — EO)M (= Re, (4.15) 
6kp CX; / 
where M3 = | wii(u) du. 


It may be noticed that the actual rate of destruction of temperature 
fluctuations by radiative transfer is 
8a TF" (0 , {08 \? 
= 2 9 ) vp a ’ 
Spc, CX; 


equivalent to an additional ‘radiation’ heat diffusivity of 


eee 

ke (4.16) 
; 3p°C 
JP p 


This was the same diffusivity as was found by Goody (1956). 


5. APPLICATION TO MOTION IN THE ATMOSPHERE 

‘This analysis is expected to apply to flows at a sufficient distance from 
solid boundaries to be substantially unaffected by their presence, a condition 
which would be satisfied in the atmosphere and outside the earth’s boundary 
laver. Within the boundary layer, the theory might be capable of repre- 
senting the order of magnitude of the effect of radiation on the turbulent 
motion, but the essential inhomogeneity of the motion would make definition 
of flow averages very uncertain. 

The logarithmic cooling rate, whose ratio to the mean velocity gradient 
determines the influence of radiative transfer on the motion, is a function 
of the scale of the temperature fluctuations as well as the radiative properties 
of the atmosphere unless the atmosphere is effectively transparent to 
temperature radiation of the gas. The condition for this is that the 


ibsorption length 
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should be large compared with the scale of the temperature fluctuations, 
which must be roughly equal to the scale of the mean velocity variation or 
the vertical extent of the flow. For a line absorption spectrum, 

2 


| 
* pk, 





where &,, is the maximum value of the absorption coefficient and may be 
very much larger than k, the ratio depending on the ratio of line spacing to 
line width. For air of sea-level density containing 2°), by volume of 
water-vapour, 


LE = 400 = cm, 


Ym 


and only motions of very small scale could be described by the ‘transparent ’ 
gmcm~%, which occurs at heights 


approximation. For air of density 10~‘ 
around 70 km, containing 2-5 x 10-4 by volume of carbon dioxide, 


v 





L = 10! km, 


> 
*m 


and a possibility of using the transparent approximation exists. 
From information given by Elsasser (1942) and by Curtis & Goody 
(1956), the mean absorption coefficient at 300° K is 


wn 


 — ( 5 15 Peel ? 
k = (90q,, + 150g, +215q¢.)gm~'cem (5.2) 
where q,, is the proportion by volume of water-vapour, g, is the proportion 
by volume of carbon dioxide, and g. is the proportion by volume of ozone. 
Using equation (4.10), the logarithmic cooling rate in transparent conditions 
at 300° K is 


Ji 
o>) 


3 = (0-21g,.+0-35¢,+0-51q.) sec}. (a 


(The rate at other temperatures can be estimated by neglecting the variation 
of k with temperature and supposing f to vary as the cube of the absolute 
temperature.) In table 1, the gradients of potential temperature, 
eT cz+g/c,, that are just sufficient to allow turbulent motion are listed for 
three values of the mean velocity gradient and for four values of the 
logarithmic cooling rate. ‘These values refer to (i) a non-radiating 
atmosphere, (ii) one containing 2-5 x 10~4 of carbon dioxide, (ii) one con- 
taining 2-5 x 10-4 of carbon dioxide and 2 x 10-3 of water-vapour, (iv) one 
containing 2:5 x 10-4 of carbon dioxide and 2 x 10~? of water-vapour, all 
effectively transparent and at 300 K. For the least absorbing atmosphere, 
the critical temperature gradient is very little different from that for a non- 
absorbing atmosphere, but for the most absorbing atmosphere, the critical 
gradient must be an order of magnitude greater. The purpose of the table 
is to show the order of magnitude of radiation effects in air of various 
compositions and not to assert that air of these compositions exists at 
great heights. The relevance (or otherwise) of these assumed compositions 
may be explained very briefly. The assumed content of carbon dioxide is 


% 
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about that known to exist in the stratosphere, and the two assumed contents 
of water-vapour have no more justification than (i) that these concentrations 
could exist without condensation at an air density of 10~’ gmcm~? and any 
likely temperature, (ii) that the occasional presence of water-vapour in 
significant quantities may be indicated by the occurrence of noctilucent 
clouds at great heights, and (iii) that these compositions lead to values 
of 8 suitably spaced for a table intended to be illustrative rather than 
definitive. 





Critical gradients of potential temperature (deg. km ') 
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) liars rai | parry 
‘ eU) - 0-005 sec} | l |= 0-Olsec! | | c {= 0-02 sec“! 
(sec”') ||ez [ez | | Jez | 
= res 2a FE |e eee 
0 | 0-064 0:255 | 1-02 
0-9 x 10-4 | 0-079 0-284 | 1-07 
5 3 e 
Stxi0* | 0-166 0-435 1:36 
43x10-3 | 1-32 | 2-63 5-26 





Table 1. 


The failure of the transparent approximation in the lower atmosphere 
makes it desirable to obtain an estimate of the logarithmic cooling rate for 
intermediate conditions between transparency and opacity. ‘The general 
expression for the factor F (equation (4.13)) may be integrated by parts 
to give 

E'(pr) dO(r) 


F , st ; 5 
| | T . ad (5.4) 





showing that F is the mean value of E’(pr)/k, taken over all values of r 
with a weighting factor dO(r)/dr. It is characteristic of line absorption 
that the initial, very rapid decrease of E’(pr) is followed by an extended 
region of very slow decrease, while the weighting factor dO(r)/dr is small 
nearr = 0. It follows that if we define a length L, by 


F = E'(pL,)/k, (5.5) 


it will be roughly equal to the scale of the temperature fluctuations. On 
this basis, critical gradients of potential temperature have been calculated 
for flow in air containing 2°/, of water-vapour at sea-level density for a 
range of scales and velocity gradients (table 2). For this atmospheric 
composition, 


PF =xFJL-' (L, in cm) (5.6) 


ustmg results quoted by Elsasser (1942). Except for the largest scale of 1 km, 
substantial effects of radiative transfer are predicted for this composition 
and mean velocity gradients of the kind that might occur at heights over 
20 metres. 
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SUMMARY 

‘The plane steady decelerated flow of a dust-gas mixture is 
analysed in an approximate manner. ‘lhe problem, which has a 
tive-parameter family of solutions, is reduced to a form such that 
the analysis can be completed by the integration of a first-order 
non-linear differential equation and a quadrature. A few integral 
curves are given and the characterizing features of the flow field 
are discussed. 


1. INTRODUCTION 

When a shock wave is propagated through a gas which contains an 
appreciable amount of dust, the thickness of the wave, the pressure change 
across the shock, and the other features of the flow differ greatly from 
those which arise when the shock passes through a dust-free gas. We 
consider here the stationary plane shock configuration which arises in a 
dusty gas and determine the dependence of the flow field on the parameters 
of interest. Since the solutions of the problems of interest are described 
by a many-parameter family of functions, no attempt is made to compile 
comprehensively the appropriate numbers. Instead, a few examples are 
worked out and the integrations which lead to detailed results are specified. 


2. "THE SHOCK PROBLEM 

Consider a homogeneous mixture of a gas and small solid particles 
which is moving with the uniform velocity uv) as in figure 1. Let the flow 
configuration be one-dimensional and let the sound speed a, of the gas 
be less than wp. It can be anticipated that a compressive change of state 
can occur without spoiling the one-dimensionality of the flow; in particular, 
when such a compressive change occurs, the first event will be a compression 
and velocity decrease in the gas consistent with the usual shock relations 
for the gas; this compression occurs in a very few mean free paths and the 
solid particles cannot be so quickly decelerated that they could modify 
these relations. Behind this conventional shock (we shall refer to it as the 
gas shock) the velocity of the gas is smaller than that of the dust and the dust 
will then be decelerated. The dust will also accept heat from the gas since 
the gas temperature has been increased above that of the ambient mixture 
by the gas shock. ‘Typical distributions of the gas velocity u(x) and the 
dust velocity 7() are shown qualitatively in figure 1. 
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| The flow configuration far downstream of the gas shock will be a steady 
one in which the gas and dust will achieve the same velocities and 
temperatures. ‘This final state can be computed very simply in the same 
way that the state following a gas shock is computed. We first write the 
equations implying the conservation of momentum and energy in the steady 
one-dimensional flow of this composite ‘fluid’. In the statements of these 
laws we ignore the partial pressure of the dust and we restrict our analysis 
to shocks which are not so intense that any evaporation of the dust takes 
place. We also require the shock intensity to be one for which the state 
law p = pRT 

is a satisfactory approximation. With the notation, p(x) = gas density, 
m = pu = const., 7(v) = dust mass per unit of volume, m = nv = const., 
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Figure 1. Qualitative gas and dust velocity distributions for the steady plane flow 
of a dust-gas mixture. 





T(x) = gas temperature, 7(x) = dust temperature (we ignore temperature 
variations within the particle), and p(x) = gas pressure, these conservation 
laws are: mu +nv +p = const., (1) 
imu? + inv? + me, T+ncer = const. (2) 
The quantities ¢,, and c¢ are, respectively, the specific heat at constant 
pressure of the gas and the specific heat of the solid material. Since pu 
is constant, the equation of state can be written 
pu = mRT. (3) 
In the initial state which we characterize by the subscript zero and in 
the final state (subscript 2), w= v and 7 = 7. Equations (1), (2) and (3) 
then imply that 
(m+n)us+mRT,/uy = (m+n)uy + mRT)/uUy (4) 


and 


4(m+n)us+ (me, +nc)T, = }(m+n)u? + (mc, +nc)To. (5) 
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Thus 2moR a. Wd. 
Uy Ug = ————- T, + —— uv, (6) 
(m+n)(o+1) a+ 
where y+nc (mc,,) 
a 
l+nc (mc) 
and y= c,/c,. However, u,, the gas speed just after the gas shock, is 
given by 
2yR y-l ‘ 
U, Uy = (a*)? = —— T, + - ue. (7) 
y+ y+1 


\ comparison of (6) and (7) shows that the tinal state and the state 
just after the gas shock will differ increasingly with increasing values of 
nc (mc,,), the ratio of the heat capacity per unit volume of the dust to that 
of the gas. 

The manner in which the state variables vary with x is governed by 
equations (1), (2) and (3), and by the specification of the mechanism whereby 
momentum and energy &re transferred from one medium to the other. 
In view of the variety of particle shapes and sizes to be anticipated in 
problems of interest, it would be optimistic indeed to attempt a description 
of these transfer processes with the size and shape dependence accounted 
for in detail. It is more profitable to write down the macroscopic rules: 


Avv, = Cpl? p(u—v)?, (5) 


Acut, = Nuk(T—r)l. (9) 


These state that the knowledge of the force on a particle of mass A and 
‘radius’ / requires only the specification of a drag coefficient C,, and that 
the corresponding knowledge of the heat transfer to such a particle requires 
the specification of a Nusselt number Vu. In general, C, and Nu depend 
on the geometry of the particles, the Reynolds number associated with the 
motion of the dust relative to that of the gas (i.e. p(u—zv)/ yu, where pu is 
the gas viscosity), and the Mach number M of that relative flow 
(\J = (u—vwv)/a, where a is the gas sound speed). However, the Mach 
number dependence is important only when is greater than unity over 
an important portion of the flow field. In view of the restrictions on the flow 
speed specified above and the increase in sound speed across the gas shock*, 
the neglect of a Mach number dependence of C,, and Nu does not imply 
any significant further restriction. Furthermore, the ratio of equations (8) 
and (9), namely 

_dr/dx 2Nuk(T—7) 
* dedx Cp lp(u—v)? 





‘ (10) 


involves only the ratio Nu/C,. It is an empirical fact that, for Reynolds 
numbers Re up to several thousand, Vu ~ KC, Re, where K is a constant 
of order unity. When this fact is used, (10) becomes 


ee 
_2 


U, u—VWV 


; (11) 





* Note that, with y = 1-4, this Mach number is only of order 2 for very strong shocks. 
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where 8 = 2Kk/(cu). In particular, if the solid material is silica, the gas 
is air and we use the values of K which are appropriate for spherical or 
short cylindrical particles, 8 = 1 is an excellent estimate. 

It is now convenient to regard v as the independent variable of our 
problem, to note that equations (1), (2) and (3) define T and 7 as given 
functions of u and vz, and to regard x(v) as a function to be determined. 
When this is done, equation (11) can be written as 





dr(u[v],v) — T[u(v), v] —7[u(v), v] (12) 
dv u(v)—v ; - 


Since T and 7 are known functions of their arguments, equation (12) 
is an ordinary differential equation from which we may determine u(zv) 
once a suitable boundary condition is prescribed. The initial condition 
which must be used requires that u(v)) = u,, the quantity defined in 
equation (7). This merely states that just behind the gas shock u = u, 
and v = uy. Equation (12) has a singularity (a focus) at the point (u,, us). 











Ue Us V 


Figure 2. The integral curves of equation (12). 


In figure 2 some of the integral curves are sketched qualitatively. That 
curve which extends from the point (,,u)) to the point (u.,u.) is the 
desired solution of equation (12). That this curve can be obtained accurately 
by any sensible numerical procedure is readily seen. ‘The nature of the 
singularity and the fact that T and 7 are quadratic in uw and @ assure a smooth 
solution curve. Once u(zv) has been determined, equation (8) may be 
integrated directly to give x(z). 

It would require an elaborate table of numbers to specify in detail the 
results of these integrations (e.g. u(x)) with regard to their detailed 
dependence on the independent parameters mn, ¢/C,,, pot @/Mor U/C, To, 
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y, and 8. Instead we shall display a few integral curves, u(zv), and record 
the dimensionless forms of the equations which would be convenient in 


computing further results. 


3. "THE DIMENSIONLESS EQUATIONS AND RESULTS 
Since no exhaustive details have been recorded, we compile in this 
section a dimensionless formulation of the problem which it is convenient 
to use for computational purposes. Let 
W = Ua; s = U/A, 6= T/T, rosea gl (Re 0=C/C,, 
€ = n'm, A = (1l+e)w,)+(yw,)7, B = (1+e)w?+(y-1)7. 
Quantities with subscript zero characterize the initial state, those with 
subscript unity the state immediately behind the gas shock, and subscript 2 
the final state. §, as defined in § 2, depends on 4, the geometric constant K, 
and the Prandtl number c, «/k. o defined as in §2 becomes 
c= (y + 5) (1 +€6),. 
With the substitution of these dimensionless quantities, equations (1), (2) 


and (3) become: 








O(w,z) = yw[A —w-es], (13) 
ae. Y ) , 
b(w, 2) ) B- — leg? — by? (14) 
€0 y= 1 _ 
L e=l,y=14,8=10, 8210 4 
: L W, 21.3, 1.6, 2.0, 3.0 
| 4 
o-—— +4 —$$$___}_—_— —_—— 
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Figure 3. Integral curves of equation (16) with « = 1. 
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Equation (12) becomes 





and, using equations (13) and (14), this reduces to 
(w—z)[yA — (y+ 1)w—yes]w’(z) = e(w—2)[yw—(y—1)-2]- 
— Be d[O(w,z)—d(w,z)]. (16) 


With the initial condition, 


W(Wy) = W, = ; 





equation (16) can be integrated to give w(z) in the interval wy > 3 > @p. 
Note that w(z) depends on the independently chosen parameters y, e, 4, 8, 
and wp. 
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Figure 4. Integral curves of equation (16) with « = 5. 
Equation (8) can now be written as 
ad 2Auy ( z dz = 
(3) = —— : (17) 





Mol J», (w—z)ReCp(Re) 
Here, A must be interpreted as the average mass per particle and, because 
of the temperature-dependent gas properties and the variable (w—z), 
Re is a function of z. For a given dust substance, A is proportional to /% 
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and for moderate Reynolds numbers C,, is proportional to (Re)-!?. Thus 
a measure of the thickness of the transition region is given by 
X = A(u>/yo 1)(uy — u,)Re'?, (18) 

where Re is evaluated at the initial state just behind the gas shock. Since X 
is just u> times the mass per particle divided by the drag per particle at the 
initial state (i.e. uy times the time required to decelarate an amount uy 
under the initial drag) the result is not surprising. Ifthe dust were composed 
of 1-micron particles of silica, for example, and the initial air flow had 
Mach number 2 and a sea-level stagnation state, XY would be about 1 cm. 
For 10-micron silica, the estimate becomes 30 cm. 

A few of the integral curves of equation (16) are drawn in figures 3 and 4 
for the values of the parameters specified in the diagrams. 
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SUMMARY 

Four alternative theoretical treatments of ‘displacement 
thickness’, and, generally, of the influence of boundary layers 
and wakes on the flow outside them, are set out, first for two- 
dimensional, and then for three-dimensional, laminar or turbulent, 
incompressible flow. They may be called the methods of ‘ flow 
reduction’, ‘equivalent sources’, ‘velocity comparison’ and 
“mean vorticity ’. 

The principal expression obtained for the displacement 
thickness 6, in three-dimensional flow may be written 


, 1 ate 
6, = 6, Uh, ay 13 
if, as orthogonal coordinates (x,y) specifying position on the 
surface, we choose x as the velocity potential of the external 
flow, and y as a coordinate, constant along the external-flow 
streamlines, such that /,, dy is the distance between (x,y) and 
(x,y +dy); and if also 6, and 6,, are the streamwise and transverse 
‘volume-fiow thicknesses ’ 

0, = 2 | (U—u) dz, 5, = é | wdaz, 
U Jy U Jy 

1s the distance from the surface, u and v are the x and y components 
of velocity, and u takes the value U just outside the boundary layer. 


1. INTRODUCTION 

A boundary layer causes the irrotational flow outside it to be that about, 
not the solid surface itself, but a surface displaced into the fluid through 
a distance 6,, the ‘displacement thickness’ of the layer, whose value at any 
point of the surface can be calculated, to a first approximation, directly 
from the velocity profile of ordinary boundary-layer theory. 

The displacement-thickness theory is well understood in two-dimensional 
flows, but it has been studied only a little in the general three-dimensional 
case (Moore 1953; Dunn & Kelly 1954). A specially useful formula for 6, 
is omitted from these papers, which select one particular method for attacking 
the problem; actually, at least four are available, all of which help to 
illuminate it, some making the theory more rigorous—less of a ‘hunch’ 
—than do most published presentations. For these reasons, the author 
felt justified in compiling a new account of the subject. 
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The discussion is limited to the case of incompressible flow, in order 
to achieve the utmost clarity of presentation. Note that, although some 
generality is lost by this, much is gained in return, since the arguments 
then apply without change both in steady and unsteady flow, and, also, 
both in laminar and turbulent flow, provided that in the latter case all words 
like ‘flow’, ‘velocity’, ‘vorticity’, etc. are taken to signify mean values, 
in the usual turbulence-theory sense. (This is because all the equations 
used are linear in those quantities.) ‘To make up in part for the restriction, 
we state without discussion, in equation (22), the form which the principal 
formula for 5, takes in the case of steady compressible laminar flow. 


‘TWO-DIMENSIONAL FLOW 
The simple case of two-dimensional flow is first used to illustrate the 
four main approaches to incompressible-flow displacement-thickness 
theory, which we may call the methods of ‘flow reduction ’,‘ equivalent 
sources ’, ‘velocity comparison’ and ‘mean vorticity’ 


2.1. Flow reduction 

With x as distance, measured along the surface, from the point of 
attachment of the boundary layer, z as distance from the surface, u and w 
as the corresponding velocities, and U as the value of wu just outside the 
boundary layer, the difference U—u represents the reduction in flow 
velocity due to the presence of rotational flow in ate boundary layer. The 


total reduction in volume flow per unit span is | (U—u) dz. 


Now, between the surface and any streamline just outside the boundary 
layer, there must be a constant volume flow per unit span. This will be 
so if the flow reduction inside the layer is compensated for by an outward 
displacement of such a streamline through a distance 6, (which produces 
a flow increase U’6,, since the velocity is U in the region of streamline 
displacement), provided that 


5, =a (U—u) dz. (1) 

This displacement of the irrotational-flow streamlines implies that they 

can be regarded as streamlines of the irrotational flow around a surface 
displaced into the fluid through a distance 6,, as stated in $1. 

We may observe that streamlines in the wake are similarly displaced, 

and behave as if the wake were a solid slab (joined to the body), of thickness 


4=s| (U—u) ds, (2) 
with irrotational flow around it. In the wake (at least in steady flow), 
5, decreases downstream at first, since the ‘momentum thickness’ 


—_— 


(U—u)u dz 


— 
w 
~- 





ui 
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must remain approximately constant; far downstream, 4,, like 4, tends 
to D/(pU*), where D is the drag per unit span and p the density. 


2.2. Equivalent sources 

The effect of the boundary layer on flow about aerofoils, particularly 
those derived by conformal mapping, can often be represented most 
conveniently, not by an effective thickening of the aerofoil, but by means 
of an equivalent surface distribution of sources. ‘lo derive this, we consider 
the form of the normal velocity zw just outside the boundary layer. In this 


region we have 


Ta on dz — | on ds ~ a ba | (U-u) dz 
oz J Ox dx Ox) 4) 
dU d [ : 
— —z+ — | l dz, + 
dx dx ( "} (4) 


since 3 is large enough for U—u to vanish, and hence for z to be replaced 
by ~ in the last integral in accordance with the conventions of boundary- 
layer theory. 

In (4), the first term is that which would be present in the irrotational 
flow around the body, and the second is an additional outflow due to the 
boundary layer. ‘This additional outflow is exactly ‘as if’ the irrotational 
flow around the body were supplemented by the effect of a surface 
distribution of sources, whose strength (volume flow rate) per unit area is 


d *% 
dx in 


—_~ 
~~ 
~ 
= 
— 
~ 
~ 
2 
\| 
| 
—_~ 
oo 
~ 
— 
\ 
J 
~— 


m 


The same considerations are applicable in the wake, where however m 
is negative, so that we may consider the rear dividing streamline of the 
irrotational flow as dotted with sinks. ‘These produce a slight acceleration 
in the mainstream flow about (for example) a flat plate at zero incidence, 
which is exactly responsible for the addition, 8-2/R, to the Blasius expression 
2-64/R!? for laminar-flow drag coefficient, which was found by Kuo (1953). 

Now, the ‘new’ fluid emitted at the sources just described would fill 
a region, adjacent to the body, of thickness 6,; for the flow of ‘new’ fluid 
past any point (with velocity U’) must equal the total outflow from the part 
of the surface between that point and the point of attachment, and this is 


[" m dx = US, (6) 


per unit span. But the external flow can be regarded as the irrotational 
flow about the surface of separation between the fluid from upstream 
and the ‘new’ fluid from the sources (as in the theory of Rankine bodies), 
and this has been shown to be a surface displaced into the fluid through a 


distance 46). 


bo 


F.M. 
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2.3. Velocity comparison 
The ‘velocity comparison’ method has some similarity to that just 
discussed, but approaches more directly the problem of finding a surface 
z = 6,(x), the irrotational flow about which is the same as the flow outside 
the given boundary layer. ‘The boundary condition for such an irrotational 
flow is 
w= l 5,(x) at z = 0,(x), (7) 
since to the required approximation the change in U may be neglected 
when multiplied by 6)(v). Hence, for values of z just outside the boundary 
layer (but small on the scale of the irrotational flow), 
w = (w) y Me 0;) U8; (x) sé 0;). (5) 
ez dx 
Comparing this with the actual form (+) of w just outside the layer, we 
obtain a differential equation 


, af . 

—(Us,)=—| (U-un) ds, (9) 
dx 

of which equation (1) represents the only solution satisfying U6, = 0 at 

the point of attachment. 


2.4. Mean vorticity 

A fourth, and perhaps most rigorous, approach regards the problem 
as that of finding the flow induced, in the presence of the body, by a given 
distribution of vorticity—or one that is, at least, approximately known. 
By ‘induced flow’ in general, one means the unique velocity field, with the 
given vorticity (and velocity at infinity), which has zero normal velocity 
at the surface; but the actual vorticity distribution is one of those, for 
which also the induced tangential velocity vanishes at the surface. 

The idea of the method is to replace the vortex layer by a vortex sheet, 


at the mean distance 





| sndz | s(cu/ez)dz | (U-—u)ds 
4 ee ra = 6, (10) 
| ds | (eu/éez) dz U 
of vorticity from the surface. Here, the boundary-layer approximation 
ou dw . cu 
= =~ ES (11) 
Cs Cox CS 


to the vorticity 7 has been used. The error in the induced flow, produced 
by such a redistribution of the vorticity at any station in the boundary layer, 
is shown below to be of the order of the square of 5// (ratio of boundary-layer 
thickness to body dimension) except near that station itself (where by 
‘near’ we mean ‘within a distance of order 6’). It follows that the 
replacement reduces the flow to rest (if terms of order (5//)? be neglected) 
at any point between the vortex sheet and the body, since the flow induced 
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at the point by ‘far’ vorticity is unaltered, and the effect of shifting all the 
‘near’ vorticity that formerly lay between the point and the surface to the 
other side of the point simply reduces the velocity there to zero. Since, 
therefore, the fluid is at rest on one side of it, the vortex sheet is a stream 
surface, and the flow outside it is the irrotational flow around that surface. 

The method gives just as directly the ‘equivalent source’ result. 
For two parallel line vortices with equal and opposite circulations K 
and — AK are equivalent to a distribution of normal doublets, of uniform 
strength A per unit area, in the strip bounded by the two lines. Hence the 
flow differs from that got by replacing all the vorticity in the layer by equal 
vorticity on the solid surface (where it has no effect, being cancelled by its 
image) by the flow induced by a volume distribution of doublets, parallel 
to the surface, whose strength per unit volume at a distance z from the 
surface is 


| nds =U-u, (12) 


since all vorticity farther from the surface than = contributes. 
To a first approximation, the effect of all these doublets is the same as 
if they lay on the surface, with strength 


| (U-u) dz (13) 


per unit area. But such a distribution of tangential doublets on the surface 
is equivalent to a distribution of sources, of strength 


@ | (U—u) dz (14) 
dx }, 
per unit area, as in (5). 

Note that the error in replacing the doublets at any station x by a doublet 
on the surface is equivalent, except very nearby, to the flow induced by a 
quadrupole of strength 


| 2(U—-u) dz, (15) 
which 1s of order (6//)*. ‘The same estimate can be similarly derived for the 
error in the original replacement of the vortex layer by a vortex sheet at 
the mean distance 6, of vorticity from the surface. 


3. "THREE-DIMENSIONAL FLOW 

Each of these approaches to displacement-thickness theory will now be 
applied in the general three-dimensional case. We take z as distance from 
the surface as before, and (x, v) as any orthogonal system of coordinates 
on the surface, such that the distance between the points (x,y) and 
(x+dx, y) 1s h, dx, while that between (x, y) and (x, y+dy) is h, dy. 
The velocities in the x,y,z directions are u,v, zw, and the values of u and w 
just outside the boundary layer are U and V. 


2B2 
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In most of the methods it is a real advantage to choose (x, y) as ‘ external- 
How coordinates ’—that is, such that the curves y = const. are streamlines 
of the external flow, which in symbols means that V = 0. ‘This implies 
that the curves x = const. are equipotentials, which suggests choosing wv as 
the velocity potential of the external flow, in which case h, U-!, However, 
we shall not require that this precise choice for x, or any particular choice 
for y, be made. 

lo fix the ideas, we assume throughout that the boundary layer attaches 
itself to the surface at a stagnation point of the external flow, from which 
all the external streamlines issue. We take x = 0, without loss of generality, 
at that point. ‘The assumption must be correct on a smooth surface, except 
in unusual circumstances like those of two-dimensional flow, when 
attachment occurs at a whole line of stagnation points (although, in two- 
dimensional theory, such lines normal to the flow plane are usually referred 
to as ‘ points’); and even then we can take x = 0 at each, since such a line 
is necessarily an equipotential. 

Note that cases when the external streamlines on the surface fall into 
two groups, those of each group issuing from ditferent stagnation points 
(such a pair of ‘nodes’ of the external streamlines being commonly 
separated by a ‘saddle-point’), are not really excluded, since the two 
parts of the surface covered by the two groups of streamlines can be treated 
separately. 

On the other hand, on surfaces with cusped edges (including the case 
of an ideally thin sharp-edged plate), when boundary-layer attachment 
occurs at the edge, it is possible for the streamline y = const. to become 
attached for a value of x, say x = x)(y), varying arbitrarily with y. We 
content ourselves with stating that all the formulae of $3 can be proved 
to be still correct in such cases if the lower limit .« = 0, wherever it occurs 


in an integral, is replaced by xo(¥). 


3.1. Flow reduction 

Consider now the portion of surface lying between two neighbouring 
external streamlines, with the constant values of y on each differing by dy. 
The boundary layer reduces the volume flow in the x-direction over this 


portion of surface by 


(h, dy) | (U-u) ds. (16) 
If the external streamlines are displaced outwards by an amount 6,, this 
makes a compensating increase (f,, dy)U6, in the said volume flow. 

However, in the general three-dimensional case, these two cannot be 
simply equated. This is because the boundary layer causes some flow at 
right angles to the external streamlines. ‘The flow across any one of them, 
between the point of attachment and the point (x, v), is 


* dx | v dz. (17) 
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This reduces the volume flow, between streamlines specified by y-coordinates 
y+dy and y, by an amount equal to the difference between the values of (17) 
on the two streamlines, namely 


(= | h,dx | @ as) dy. (18) 


ey “0 “UU 





Accordingly, part of the total flow reduction (16) is associated with this 
cross-flow effect, rather than with displacement of the external streamlines, 
whence (f, dy)U6, must be the difference of (16) and (18), giving 


; Lf : 1 oa ff i 
6,=— | (U-u)dz— —-—| h,dx| vwds. (19) 
U }, Uh, cy J Jo 
In terms of the streamwise and transverse ‘ volume-flow thicknesses ’, 
defined as 
Lot | Gwe, ual om 20 
0, = U } ( —u) ds, 0, = U | C 43, (2 ) 
we have 
| es ay ee 
(a) = 0.- —=—_- = Uh ..6, dx. (21) 
. Uh, ey | - . 


y & 
If x is taken as the velocity potential of the external flow, the product Uh, 
in (21) is 1, which further simplifies the equation. 

The same method shows that the expression analogous to (21) in the 
case of steady compressible laminar flow is 


— BUR. = | PUh,8, dx, (22) 
where P is the value of the density p just outside the boundary layer, and 
6, 6, are the ‘mass-flow’ thicknesses 

6. = es | (PU-—pu) ds, $= ea | pv dz. (23) 
JO Je r PUI, 


We shall not, however, consider compressible flow any further. 


3.2. Equivalent sources 
For this method we need the equation of continuity in boundary-layer 
coordinates, namely 
1 7 C Cw > 
——< — (h, u)+ —(h,v)>+ — =9, (24) 
hh, \ ex ey dz 


4 


from which it follows that, just outside the boundary layer, 


Ow ‘* — ie 0 


—(h,u)+ —(h,v)$ ds 
9 





- J 


= - FE lO) mon { ~u) ds) x(k. { eds)}, 
(25) 


where = has been replaced by ~ in the integrals because it is supposed large 








enough for U—u and ¢@ to vanish. 
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In (25), the first term is that which alone would be present in the 
irrotational flow around the body. ‘The rest is the additional outflow due 
to the boundary layer, and is ‘as if’ there were a source distribution on 
the surface, of strength 


~~ (0%. 8,) — (Uh,4,)> (26) 
Ox 


hh, 


per unit area. In many applications, one could usefully consider the 


m 


( 'y 


irrotational flow around the original solid surface as modified by this 
‘equivalent source’ distribution. Alternatively, if we wish, we can deduce 
the formula for 6, from it, as in $2.2. For ‘new’ fluid is created at these 
sources, in the region between streamlines of the irrotational flow with 
y-coordinates y and y+ dy, at a rate 


( m(h, dy)(h,.dx) = dy{ Uh, 6,.- = | Uh,.6, dx }. (27) 
"ea 4 / Oy / 


Hence, if 5, is the thickness of the layer of ‘new’ fluid at the point (x, y), 
then (4, dy)U6, must equal (27), which gives equation (21) for 6). 


3.3. Velocity comparison 

The ‘velocity comparison’ method has already been described by 
Moore (1953), in general orthogonal coordinates (x,y). We therefore 
confine ourselves to an account of the method in the special case when 
external-flow coordinates are used. 

If s = 6,(x, y) is a surface, the irrotational flow about which is the same 
as the flow outside the given boundary layer, then in that irrotational flow 
the boundary condition is 

Ww ad — at z = 0,(x, y). (28) 
Rex c 
Hence, for values of z just outside the boundary layer (but small on the scale 
of the irrotational flow), 
site ele mw) «, Ue, s=0, @ ‘i 
w= (w), 44 (3) a ee 
Comparing (29) with the actual form (25) of w just outside the boundary 
layer, we obtain (after multiplication by h,,/,,) a differential equation 
cy 





~ 
rai 


(Uh, 5,) = —(Uh,8,)— —(Uh,8,), (30) 

Ox ex 

of which equation (21) is the solution such that UA,,6, = 0 at x = 0. 
Moore’s application of this method in general surface coordinates (x, y) 

vields the partial differential equation (written by him in vector form, 

but here translated into scalars) 


=| hy Us, (U—-u) dz | } =| (V8, | (V —v) dz | = (), (SP) 
cx CY 0 


but he does not solve it for 5,. Since the external streamlines are the 


characteristics of equation (31), its solution in general must necessarily 
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involve integrating along those streamlines, as was done above. On the 
other hand, there may be particular cases when the geometrically simplest 
choice of coordinates (x, y), rather than external-flow coordinates, happens 
to give (31) a form whose solution can be spotted without any integration, 
and in such a case Moore’s approach would obviously be better. 


3.4. Mean vorticity 
Finally, we apply the ‘mean vorticity’ method to three-dimensional 
How. On the boundary-layer approximation, the v and y components of 
vorticity are 
. ev Cu . 
Ket, weet (32) 


CS C2 





Integrated across the layer, in external-flow coordinates, these give 
é dz = 0, | 9 dz = U. (33) 
Thus, the mean vortex lines (averaged across the boundary layer) lie along 
the equipotentials 1 = const. of the external flow. 

If this y-vorticity is now redistributed, by concentrating it in a vortex 
sheet at the mean distance 


sn dz | (U-u) dz 


/ 





| dz U 
of y-vorticity from the surface, the error is of order (4//)? as in § 2.4, so that 
to the order of approximation required the effect on the external flow is 
unchanged. If there were no x-vorticity we could then argue, as before, 
that we are left simply with the irrotational flow around a surface displaced 
into the fluid through a distance 64... 

However, we cannot neglect the streamwise vorticity simply because 
its average across the layer is zero. ‘This means, rather, that it can be 
regarded as made up of small vortex rings, in planes perpendicular to the 
surface equipotentials. From the equivalence of a vortex ring of circulation 
K to a shell of normal doublets, of strength A per unit area, whose boundary 
is the ring, it then follows, as in §2.4, that the streamwise vorticity € is 
equivalent to a distribution of doublets with axes in the negative y-direction, 
of strength 


Edz=v (35) 


per unit volume. ‘These doublets are equivalent, with an error of order 
(6 /)?, to a surface distribution of doublets of strength 
/ 


| wvdzs= Ub, (36) 


per unit area. ‘This distribution of tangential doublets with axes in the 
negative y-direction (geometrically, around equipotentials) is equivalent 
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to a distribution of sources of strength 
1 d,,, - 
~ Fh, By (U8v'e) (37) 
per unit area. 

The induced flow consists, therefore, of the irrotational flow around 
the vortex-sheet, at distance 5, from the surface, supplemented by the 
effect of the sources (37), which may be shown by the methods of § 3.2 
to alter the effective displacement thickness from 6, to 

2 sae 
6, = 6,- —-—] Uh,6, dx, (38) 
Uh, cv Jo 
in agreement with (21). 

Alternatively, we can deduce from this approach the equivalent source 
distribution (26). For the flow differs from that got by replacing all the 
y-vorticity by equal vorticity at the surface by the velocity field of a 
distribution of doublets with axes in the positive x-direction, of strength 


( ndz= U-u (39) 


per unit volume. ‘This is equivalent to a surface distribution of doublets 


of strength 


| (U—-u) dz = U6, 


per unit area, and this distribution of tangential doublets, with axes in the 
positive x-direction, is equivalent to a distribution of sources of strength 
we 
——— — (U/5.h,) (40) 
hh, ex 
per unit area, which together with (37) gives the original result (26). 

We may conclude with a remark about the influence of wakes in general, 
three-dimensional flows, in which (as well as in unsteady two-dimensional 
flows) the vorticity integrated across the wake may have a non-zero value 
namely, the trailing vorticity. ‘The method of this section gives that the 
error in the usual theory, in which this vorticity is regarded as concentrated 
into a vortex sheet (which we may take as z = 0, with x,y as orthogonal 
coordinates specifying position on the sheet), is equivalent to the flow 
induced by a surface distribution of sources on z = 0, whose strength 
per unit area is obtained by adding the source strength (26) for the 
‘boundary layer’ on one side of the surface z = 0 to the value of (26) for 


the ‘boundary layer’ on the other side. 
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SUMMARY 


A study is made of the propagation of sound in both a constant 
gradient shear flow and a turbulent shear flow above a flat surface. 
Curves are presented showing how, in the case of downstream 
propagation, the flow gradient tends to channel the sound energy 
into a narrow layer next to the wall. ‘These results are used in 
estimating the effect of a flow on the attenuation of sound in a duct 
with absorbing side walls. 


1. INTRODUCTION 

Undesired sound is frequently attenuated by having it pass through 
ducts with absorbing side walls. If the absorptive properties of the walls 
are described by a normal impedance, and if there is no air flow in the duct, 
the resulting attenuation is easily predicted on the basis of existing theory 
(see, for example, Morse 1939). In many cases, however (e.g. in ventilating 
ducts or the exhaust ducts in certain wind tunnels or jet-engine test cells), 
there is also a flow of air through the duct which may affect the rate of 
sound attenuation. Downstream, for instance, the velocity gradients in 
the boundary layers might be expected to direct the sound into the absorbing 
walls and so increase the attenuation. 

A particularly simple case to study is that of sound propagating in a 
fluid which is flowing between two parallel walls. Here the flow velocity 
can be assumed to be approximately constant throughout a central region 
between the walls and to fall to zero within the two boundary layers. Since 
the flow profile will be symmetric about a centre line drawn midway between 
the walls, it is then sufficient to consider just half the region, i.e. that between 
the centre line and one wall. On a ray acoustics picture it is clear that 
sound rays propagating upstream in one of the boundary layers will be bent 
away from the wall, while rays propagating downstream will be bent towards 
it, and in certain cases may be expected to suffer repeated successive reflections 
in the wall surface. A ray treatment is not, however, generally adequate for 
ascertaining the actual sound pressure distribution across such a shear layer 
owing to the tendency, pronounced at lower frequencies, for the acoustic 
energy to ditfuse away from regions where it is concentrated. It is therefore 
desirable to adopt a wave treatment from the outset. 

* Now at Mechanical Engineering Department, Massachusetts Institute of Tech- 
nology 
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We shall confine our attention to the problem of a two-dimensional 
1coustic wave propagating downstream in a shear layer above a plane wall. 
Within the shear layer we shall consider only the lowest mode of propagation. 
Thus the boundary condition to be imposed on the sound wave will be that 
the normal component of the particle velocity (in a direction perpendicular 
to the wall) must vanish both at the wall surface and at the top of the shear 
layer, so that the wave matches smoothly on to a plane wave above the layer. 
First we study the simplest case of perfectly rigid side walls and a constant 
gradient of flow velocity, and then we extend the analysis to a turbulent 
How profile in which the flow velocity increases as the one-seventh power 
of the distance from the wall. For both these cases representative curves 
are given showing the variation of sound pressure across the duct which is 
brought about by the presence of the flow gradient. Finally, we shall 
study the effect of the shear flow on the attenuation of the sound when the 
walls have a small admittance. In this case the flow has two contrary 
etfects on the sound transmission. In the first place (for downstream 
propagation) it tends to increase the absorption by directing sound into 
the walls as previously pointed out. On the other hand, for wavelengths 
long compared to the shear layer thickness, this refraction by the flow 
gradient becomes quite negligible and may in fact be counterbalanced by 
the increase of intensity in the central region which results from the 
transport of acoustic energy by the flow. 


2. THE BASIC EQUATIONS 
We shall first derive the wave equation in a form suitable for studying 
the propagation of sound in a shear flow, If we consider the flow velocity 
V’ to lie in the x-direction and to be a function of y only, then, neglecting 
viscosity, we obtain the linearized Navier-Stokes equations for the 
conservation of mass and momentum in the form 


t OX ox C 
+i a0: co a (2) 
ct CX O71 Py CX 
ae oe (3) 
ct Xv Po CY 


Here py represents the static density of the medium, which is assumed 
constant; the density fluctuations p, the sound pressure p, and the particle 
velocity components u and wv are small quantities of the first order. If the 
further assumption is made that p = cp, that is, that the pressure and density 
are adiabatically related, then these equations yield 


where VW = Vc. 
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We shall be interested in solutions to equation (4) of the form 


p = eik(xx-a) F(x, y), | 


, (5) 
U = eM«r-dG(K, y), 
where F and G are related by equation (3), 
- = Ipyckh(1—«M)G. (6) 
dy 


Substituting these expressions into the differential equation (+), we obtain 


@F  2xM dF . 
+ © +[(1—-«M)—2]F = 0, 
eae eee “) 


where M'’=dM/dy. ‘This equation must be satisfied by the pressure 
amplitude F subject to the boundary conditions dF/dy = 0 at y = 0 and 
y= L. 

Although the approximate method of solution we shall employ can be 
applied to a fairly wide class of velocity profiles, the analysis is considerably 
simplified by taking the flow to have a constant velocity gradient. Actually, 
the velocity of a turbulent fluid flowing over a flat plate increases approxi- 
mately as the one-seventh power of the distance from the plate. We shall, 
however, first consider the simple case of a constant gradient before treating 
the physically more realistic situation described by the one-seventh power 
law. 





3. 'THE CONSTANT GRADIENT BOUNDARY LAYER 
For the constant gradient case we put WM’ = k/A where A is a constant 
and introduce 7 = «~!— M as the independent variable in (7), which then 
reduces to 


@F 2dF , 





Ca. A22(n2—1)F = 0. 8 
dyn® dn i ill ) 


Approximate solutions to equation (8), valid asymptotically for large A, 


can be obtained by a method proposed by Langer (1937). ‘These are 


F = ns"*q-"* f(u), (9) 


W here 


Bs 9 2 2 . . 
s= | gq _ dy = Nn — 1)! *— cosh NS» 


u= (3Aks)? 7 
and f(w) is the general solution of 
f’(u) +uf(u) = 0. (9a) 


These solutions will, of course, break down at y = 0, i.e. when the flow 
velocity equals the propagation velocity, and hence cannot be used through 
this point. We shall accordingly assume always that the flow is subsonic 
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and 7 > 0. The function f can be expressed either in terms of one-third 
order Bessel functions, namely u!?J,3(3u3?) and u'*J_,4(2u%?), or as Airy 
integrals, 4i(—u) and Bi(—u). 

If we denote two independent solutions of (8) by F, and F;, then the 
combination of these which satisfies the boundary conditions is 


F = F,-RF,, 





where 
LF, /dy 
R= ie 
(FF Fe sounder 
and 
R(y = 0) = R(y = L). (10) 
In terms of Airy integrals we have 
F = ns! *q-'4(Bi(—u)— R Ai(—u)], (11) 
and 
, Bas -)2/3 Ry ( 
R ¢ Bi( —u) — (Ax)? Br'( od (12) 
fh Ai( —u) — (Ax)?3 Ar’(—u) 
where 


fh = (12s)-?3[1 —3(2 — y?)/(yq*?)]. 

Equation (10) determines the eigenvalues of x. We shall be interested 
in the smallest one of these, corresponding to the lowest mode of propagation 
in the duct. For particular values of A and kL the equation can be solved 
for x by a method of successive approximation. It is clear on physical 
grounds that for downstream propagation the propagation velocity cx~! 
must lie between ¢ and c(1+M,), where M, is the value of M at y = L. 
When there is little refraction by the flow, i.e. for low frequencies or small 
flow velocities, cx! will lie nearly midway between these values. On the 
other hand, in cases where the refraction is large and most of the sound 1s 
propagating near the walls, then cx~! will be nearer the smaller of these 
values. ‘These considerations are a help in guessing a first value for «. 

For purposes of calculation the quantities s, u, @ occurring in (11) 
and (12) are most conveniently expressed as power series in the small 
quantity « = 7-1 = «~!—M-1, namely, 


v2 


? 
oS" 63211 +.0-150¢ —0-134e2 +...], 





> 
Pi 


u = 213( Ax)? 3 e[1 +0-100¢ — 0-012e? + ...], 
h = 0:7143 —0-8938e + 0-955e? 4+ ..., 
F = 3-18(1+0-90e)[Bi(—u) — R Ai(—u)]. 


Figures 1 and 2 are plots of the sound pressure level in decibels (i.e. of 
20) logy)(P/p,)) across one of the boundary layers in the duct as a function 
of the distance from the wall for different values of .W,, the Mach number 
at the top of the layer. In figure 1, kL = 27, that is, the boundary layer 
thickness L is equal to a wavelength, while in figure 2, RL = 20. The 
shear flow is seen to have a much greater effect on the high frequencies 
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Figure 1. Sound pressure level in decibels (20 log, (p/p,)) as a function of distance 
from the wall for a sound wave of frequency parameter kL 27 which 


propagates in a flow whose Mach number increases from 0 at the wall to .V/, at 
1 distance L from the wall. 
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Figure 2. Sound pressure level in decibels as a function of distance from the wall 
when kL = 20. 
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than on the low; in fact, for the higher frequency (kL = 20) the theory 
predicts a pressure difference of as much as 90 decibels for a flow of Mach 
number 0-5. ‘The Mach number dependence of the sound pressure profile 
is illustrated in figure 3, where the sound pressure difference in decibels 
across the shear layer has been plotted as a function of Mach number for 














the two cases RL 27 and kL 20. 
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Figure 3. Sound pressure level difference in decibels across the shear layers of 
, figures 1 and 2. 


4. ‘THE TURBULENT BOUNDARY LAYER 
We now consider a typical turbulent boundary layer in which the average 
How velocity varies as the one-seventh power of the distance from the wall. 
Setting VW = M,(y/L)!* in (7) and using M as the independent variable 
we obtain 
d?F 6 2\ dF } - et ba td : ‘ 
ll — — -)—_ 4+40(kL)?M> 42M }2(n2-1)F =0, — (13) 
dv? \M 7») dM 
where, as before, 7) Pie eee | 
Langer’s method can again be employed to give asymptotic solutions 
to this equation in the form 
F = M*ys'®q-!4[Bi(—u) — R Ai(—u)], (14) 
where 


g = M*(7?-1), u = ($Axs)**, 


$= [ q' 2 dn, A 


l 


7RLM,’. 


ll 


The boundary conditions are satisfied, as before, by requiring that R, be 
equal to R,, the subscripts referring to values at y = Oandy = L, respectively. 
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A difficulty now presents itself in determining the form of R at y = 0. 
This arises because the coefficient of dF/dM in equation (13) is singular 
at y = 0 (i.e. .W = 0), and hence the equation does not satisfy the conditions 
of Langer’s theory right at this point. Langer’s solutions are, however, 
valid asymptotically (\—- 2%) in the neighbourhood of y=0. Now 
equation (13) has only a regular singular point at .W7 = 0, and hence its 
independent solutions can be found as Laurent series with indices 0 and 7. 
The first few terms of the general solution turn out to be 


F = ¢,(1—AA2(1—«?)M™ +...) 4 


dx 


+ ¢o(M7(1 — fe M + 2x2 M2) — 1 A2(1 — 2) M24...), (15) 


where c, and c, are constants. Remembering that y is proportional to MW", 
we see that the condition dF dy = 0 at y = 0 rules out the solution which 
starts as J. Hence we can replace the boundary condition at the wall by 
the condition that as \ — « and WM 0, equation (14) must tend to 
equation (15) with c,=0. The calculation involves substituting the 
asymptotic forms of the Airy integrals into equation (14) and expanding the 
other quantities in powers of WZ. The result turns out to be that (15), with 
Cy. = 0, is asymptotically equivalent to 
MB ys! ®q-**[ Bi( —u) + (Al( —uy)/ Bi( —u,))Al( —u)], 

where uw, is u at y = 0. Thus we see that 


R, = — At(—u,)/Bi(—w,). (16) 


At vy = L equation (13) is well behaved and hence R, has the form given in 
equation (12), where now 


Again, for calculation it is useful to expand the various quantities 
appearing in (14) as series in e = 7—1 and B=e/(x-!—1). First, 


s= 2 2(«-!— 1)%e3?(a,+a,€+a,€7+...), 
where 
Gy = 3(4 — $8 + PB? — 0B + 1884 — $85 + 2 88), 
a, = 3(3 -$8 + BB... + 48), 
a, = —3(i—38 + 2B? —... + £89), 


and so on. Secondly, 
“i= ?1 3(AK)* 3(k — 1 )a? 3(] _ Cy e+ Co € i a ace hy 


where 


| 
wil bo 
{A 
a 
! 
| 
Wo! bo 
A138 
=) 
— te 
Ol = 
— 
c“ a“ 
=> |S 
al 
t 


Thirdly, 
hb = 2-*3(n-1- 1)a~*3e-1(d,+d,e+d,e7 +...), 
where 
d, = 1-5, d, = 0:25 —3-5b, —b, —c,(1—4,), 
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and h 


Finally, 
F = 3-!8(«-!— 1)a'8(1+9,€+9,e7+...){/Bi(—u)—RAi(—u)}, 


Sl 52 
vhere 
= ae 
2, Oss -+ ()-167 il 
Ay 
a ae = fain? 
a, (0-086 + 0-146 —! + 0-167 —2 —0-07 i= 
~ Ay ay ay 


\s an aid in calculation the quantities a= *, a/®, —c, and g, have been 
plotted in figure 4 as functions of 8 = 1— M(y)/(«7'—1). Note that in 
order to evaluate F for a given set of values of the parameters .V, and kL, 
































Figure +. A plot of various quantities appearing in the numerical evaluation of 
equation (14). 

it is necessary first of all to ascertain the value of «~! which satisfies the 

boundary conditions, R, = R,, where these quantities are given in 

equations (16) and (12), respectively. In these expressions uy means the 


value of u at y = 0, ie. at e = «-!—1, and u, is the value of u at v= L 
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or » =x !'—1—M,. Having obtained «~! by a method of successive 
approximation, one then has f, and hence the various coefficients in the 
above power expansions, as functions of y, and, in terms of these, the 
function F is determined. 

Figure 5 is a plot of the turbulent boundary layer profile M/M, = (y/L)!? 
as a function of y/Z. In figure 6 the sound pressure variation across such 















































1-0 T T T T T 
== 
— ——" | “a 
08 at ns a 
06 es eae | 
M | Y 
M, | | 
04 t t 
3 | | 
02 =F t a 
| | 
a Z| 
| 
! a a | heoaial 
¢) O-2 0,4 0,6 0,8 1,0 
Y/L 
Figure 5. Turbulent boundary layer profile M/M, = (9/L)"/7 as a function of y/L. 
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Figure 6. Sound pressure variation across the flow profile of figure 5 for M, = 0-2 
and kL = 20. 
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a profile is plotted for the case M, = 0-2, RL = 20. ‘The dip in the curve | 
near y = 0 does not have a physical basis but is merely a manifestation of 


the fact that the Langer solution breaks down at this point. ‘The curve is 
qualitatively similar to that for /, = 0-05 in figure 2, which was drawn for 
the case of a constant flow gradient. ‘This is physically not surprising if 


we consider the shape of the turbulent boundary-layer profile. Referring 


to figure 5 we see that the one-seventh-power curve can be replaced with 
a fair degree of approximation over most of its range by a straight line of | 
slope /, 41 starting from M 3M, at y = 0, as is indicated by the dotted 


line. ‘The area under the two curves is the same, although the dotted line 

involves a change in M of only |.M,. ‘Thus, we see that a reasonable estimate 
of the sound pressure distribution across a turbulent shear layer is obtained 
by considering the distribution across a suitable constant flow gradient, 


for which the calculation is far less laborious. 


5. EFFECT OF FLOW ON THE ATTENUATION OF SOUND IN A DUCT 


We now suppose the side walls to be absorbing, so that the sound is 
attenuated as it propagates, and we ask what effect the presence of an air 
How will have on this attenuation. We assume that the absorption of 
acoustic energy by the walls is sufficiently small so that the sound pressure 
profile across the duct is not sensibly different from that which has already 
been calculated on the assumption that the walls were non-absorbing. 


Then the power transmitted down a unit width of the duct in one of the 
boundary layers is 
“L 
a | I(y) dy, 
0 
where /(y) is the intensity of the sound field in the absence of absorption. 
yower absorbed in an interval dv of the side wall 





On the other hand, the | 
will be approximately 


dil {Es} A(Z) dx 
_ Po x dx 
Po¢ 


Here py is the value of p at the wall y = 0; #(Z) denotes the real part of 
the wall impedance 7, and x pyc = -A(Z)) Z? is the real part of the wall 
admittance. ‘Thus the power in the duct diminishes according to the 
relation 
2 OX a 

ae (17) 
Po ¢ I(y) dy 

If there were no flow and the acoustic pressure were uniform across 
the face of the duct, then we should have / = p?/(pgc) and 

W = W,e-27Z, 


HW’ = Wyexp - 


This suggests writing (17) in the form 


W = W, e274, 





—— 
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where 


“L 7 (18) 


The quantity 7 = n(kL, M), which is a function of both the frequency and 
the flow velocity, is then a measure of the increase in attenuation due to the 
presence of the flow. In other words, if the acoustic intensity in a certain 
frequency component is attenuated at the rate of, say, 6 decibels in a unit 
distance without flow, then in the presence of the flow it should be attenuated 
by 6n decibels in the same distance. 

In order to evaluate m we must first obtain an expression for the intensity 
I(y). It is more difficult to determine the acoustic intensity in the shear 
flow than it is to obtain the pressure profile across it owing to the fact that, 
unlike the pressure, the intensity is a second-order quantity (involving the 
square of the acoustic pressure). Hence, in deriving the intensity it would 
clearly be unjustified to neglect second-order terms in the basic equations 
as was done in linearizing equations (1) to (3). 

In this paper we shall confine ourselves to deriving an expression for the 
acoustic intensity, correct to second order, for the one-dimensional case 
of a plane sound wave propagating in a medium which is moving with the 
(constant) velocity 1’ = cM in the direction of propagation. This is done 
in the Appendix, and the result turns out to be 


r= 2 (114M). (19) 


We shall then make the assumption that this expression represents a 
reasonable approximation to the actual intensity /() in the shear flow if p 
and VW are replaced by p(v) and M(y), i.e. by the known values of these 
quantities in the shear flow. This approximation is presumbly better the 
higher the frequency. 

On introducing /(\) from (19) into (18), we obtain 


“f 2 ; 
gl (4) (1-4 wy ® 
0 \Po E 


ial £ Dg , dy 
- sM)%. 20 
| (=) a+nF (20) 


In accordance with the discussion of the one-seventh-power profile 
of figure 5 given at the end of the last section, we suppose the Mach number of 
the flow to start from a value M, = 3M, at y = 0 and to increase uniformly 
to M, at y= L. Now the pressure variation across the shear layer, given 
by F, depends on the velocity gradient only and is not changed by the 
superposition of the constant flow /,. Thus the function F is still given 
by equation (9). On the other hand, the factor (1 +.) in equation (20) 
clearly involves the total Mach number. ‘Thus, to the first order ine = 7—1 
we have 


F = 3-16(1 +0-90e)f(u), 
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where f(u) = Bi(—u)—R Ai(—4), 
and u = 2'3(Ax)*3e(1 + 0-le). 
Also 1/3 
: a (=) (1—«/5) du, 


M = M,+x«1—-7. 


Substituting these expressions into equation (20) gives 








n—* = 3(«-!-1) l WF, " (1+-[u) du 
where 8 | ; 
ee cca 7-1/3 2/3 
(; i+ a) oii 04 Ali 
4 T a | | 
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Figure 7. 


0-4 at a distance 1 from the wall. 


All quantities appearing in equation (21) are obtained from the results 
of the first section applied to a Mach number $M. The integral in (21) 


can be evaluated with the help of the following relations 


| f2(u) du = uf? +f", 


| uf?(u) du = \(uf’? —ff' + uf?), 


which can be established directly from the equation (9 a) satisfied by f. 

In figure 7 we have plotted m vs RL for a Mach number of 0-4. 
clear from the figure that the effect of the flow is to increase the attenuation 
of the higher frequencies (RL > 5) but to diminish that of the lower 
In figure 8 we have plotted m vs the Mach 


frequencies (RL < 5). 


number for kL = 27 and kL = 20. Here we see that for the higher 





\ plot of the attenuation parameter » as a function of RL for M, = 0-4. 
The flow velocity increases linearly from a Mach number 0:3 at the wall to 


It is 


—— 


 ——— 


—— 
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frequency increases with Mach number up to a certain point but then 
begins to diminish again when the Mach number becomes too large. On 
the other hand, for the longer wavelength (kL = 27) n remains always 
near to unity, indicating that for this case the flow practically annuls the 
attenuation. As pointed out in the Introduction, these results are due to the 
opposing effects of refraction by the flow gradient and increased intensity 
due to the mean flow. The presence of these two factors is clear from the 
expression (20) for » where the term (/'/F,)* represents the contribution 
from refraction and the term 1+ © that from the mean flow. 





4 T T T T T 
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= 
| lL | | 
O 0:4 0:8 1-2 6 2:0 
M, 
Figure 8. A plot of » vs M, for kL = 27 and kL = 20. ‘The flow increases linearly 


from {.V/, at the wall to .W/, at a distance L from the wall. 


It should again be pointed out that the analysis in §3 has been applied 
only to subsonic Mach numbers. ‘Thus, for example, the point M, = 1 
on the abscissa in figure 8 is to be interpreted as meaning a constant flow 
with VM, = 0-75 superposed on a shear flow which increases linearly from 
M=0 to M=0-25; §3 has been required only in dealing with the 
superposed shear flow. 


The author is indebted to Professor Lighthill for several helpful 
discussions of this problem. 
APPENDIX 
Acoustic intensity in a moving medium 
The energy flow in a moving medium is given by (cf. Blokhintsev 1956, 
p. +) 

N = (!pv? + pE+p)v, (1 A) 
where p is the density, FE the internal energy per unit mass, p the pressure, 
and v the total velocity. For an ideal gas pE = p/(y— 1) where y is the ratio 
of the specific heats. We now write 

P=PotPitP» P = Pot Pit Pas v= Viv,+v3, 
where fp, is the undisturbed pressure, p, the first-order acoustic pressure, 
and p, the second-order acoustic pressure (of the order of pj) and so on. 
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If we define the acoustic intensity J as the difference between the flow 
vectors N as calculated from (1 A) with and without the sound field present, ? 
we obtain with an accuracy up to terms of second order 


[ = N—N, = §po02 V +3 py V¥0,+3p,0, V2+ 4p, V4 
tiv ly 1)} (po V+p,%+Ppov2), (2A) 
where the bar over .V— N, indicates a time average, but bars are omitted 
from terms on the right for convenience. In deriving the above expression | 
it has been assumed that the velocities V, v, and v, all lie in the same direction 
and that the time average of all first-order quantities is zero. 
We next express all second-order quantities in terms of first-order 
quantities as follows. ‘The equation of conservation of mass in one 





dimension 2 R 
co Oo 
= + = (pv) = 0 
ot Ox 
becomes on time-averaging 
Pot2+pi%,+pol = 0. (3 A) 
Similarly the equation of conservation of momentum 
0(pv) rai 5 op 
——$ ++ — (pz a) = sabe 
ot ox ox 
g1VeS Po U7 + 2p Vv, . 2p, 0; V 4 Pes V? 4 pe ), (4 A) 
Finally, the equation of state, p = Po(p/py)”, gives 
pi = "pi, (5 A) 
and ‘ y-1 45 
Po = Cp,+ >— "pj, (6A) 
where c* YPo/Po- “Po 
The three equations (3 A), (4A) and (6A) determine p., p. and v, in 


terms of zero- and first-order quantities ; 


’ 


po = — Py CF “37a aw? 


2(1— M2) 
vi Coytl 
Po = —Po 2 2(1- MP) 
ee. ae he ee | 
7 Po 2 1-M?} 


where M = Vc. 
If we now substitute these expressions into (2 A) and make use of (5 A), 
together with the plane wave relation p, = py czv,, we obtain, after reduction, 
r= Pi 44M). (7A) 

Po¢ , 
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SUMMARY 

The theory of an ‘ideal dissociating ’ gas developed by Lighthill 
(1957) for conditions of thermodynamic equilibrium ts extended 
to non-equilibrium conditions by postulating a simple rate equation 
for the dissociation process (including the effects of recombination). 
This equation contains the ‘equilibrium’ parameters of the 
Lighthill theory plus a further *‘ non-equilibrium’ parameter which 
determines the time scale of the dissociation phenomena. 

The behaviour of this gas is investigated in flow through a 
strong normal shock wave and past a blutf body. ‘The assumption 
is made that the gas receives complete excitation of its rotational 
and vibrational degrees of freedom in an infinitesimally thin 
region according to the familiar Rankine-Hugé@hiot shock wave 
relations before dissociation begins. ‘lhe variation of the relevant 
thermodynamic variables downstream of this region is then 
computed in a few particular cases. ‘The method used in the latter 
case is an extension of the ‘Newtonian’ theory of hypersonic 
inviscid flow. In particular, the case of a sphere is treated in 
some detail. The variation of the shock shape and the ‘ stand-off’ 
distance with the coefficient .\, which is the ratio of the sphere 
diameter to the length scale of the dissociation process, 1s exhibited 
for conditions extending from completely undissociated flow to 
dissociated flow in thermal equilibrium. Results would indicate 
that significant and observable changes from the undissociated 
values occur, although values for the non-equilibrium parameter 


are not, at present, available. 


1. INTRODUCTION 

Lighthill (1957) has considered the fluid dynamical problems associated 
with a dissociating gas under conditions of thermodynamic equilibrium, 
while referring to papers in preparation on the quasi-equilibrium and 
non-equilibrium conditions. In order to make the problem tractable an 
idealized gas was postulated which exhibits, during dissociation, the main 
features of a real dissociating gas such as nitrogen or oxygen. ‘This gas 
was called ‘an ideal dissociating gas’. In this paper, two problems will 
be considered in which the behaviour of this gas is studied under non- 
equilibrium thermodynamic conditions. ‘Thus, we attack the problem 
originally proposed (Lighthill 1957) for part II] of the series of which that 
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paper was part I. To do this, it is necessary to postulate the rate of 
dissociation of the molecules into atoms, under conditions tar removed 
from equilibrium. At present the form of such a relation and the values 
of the relevant parameters concerned with its non-equilibrium behaviour 
are not known, although attempts have been made to deduce such a relation 
(see, for example, Wood 1956 and Evans 1956). In view of this uncertainty 
a somewhat simplified form of the dissociation rate has been assumed 
which gives the correct equilibrium conditions, as well as the most important 
features of the non-equilibrium ones. ‘This relation is discussed in §2. 
It will be seen in the later sections of this paper, however, that the exact 
form of this relation may not give results which ditfer very much from the 
present one, owing to the predominance of one factor over all others. ‘This 
is the exponential factor in the Maxwellian distribution of the energy in 
the molecules and atoms. 

The two problems here considered are, first, the flow of an ideal 
dissociating gas, assumed initially undissociated, through a strong normal 
shock wave and, secondly, the flow of the same gas past a blutt body. The 
latter problem involves a further generalization of the author’s theory of 
How past blutt bodies on the so-called ‘ Newtonian’ assumptions. 

‘The gas is assumed in both cases to have its rotational and vibrational 
modes fully excited by passing through an infinitesimally thin shock wave 
behind which the gas begins to dissociate. In other words we assume 
the time scale for the changes in the rotational and vibrational energy of the 
molecules is small compared with that for changes in the amount of 
dissociation in the gas. When undissociated, the Lighthill ideal dissociating 
gas behaves like a perfect diatomic gas with constant specific heats but 
with only half the vibrational energy as given by the principle of equi- 
partition of energy. ‘Thus the ratio of specific heats when undissociated 
is 4. If initially, the flow has a velocity U and density py), passage through 
a normal shock wave, in the absence of dissociation, reduces the velocity 
to $l , and increases the density and pressure to 7p, and § Sp) l® respectively. 
The effect of dissociation is then to absorb some of the energy associated 
with the translational, rotational and vibrational modes of the molecule, 
thereby reducing the temperature of the gas. The pressure itself can change 
very little since the maximum it can reach is py U*. Hence, the main effect 
of dissociation is to further increase the density behind a normal shock 
wave. In the above argument ‘temperature’ is used to mean that 
associated with the energy in the translational modes of the molecules 
and atoms. 

On the assumption that the normal shock is strong, the main parameter 
that describes the flow through the initial shock wave is the energy per unit 
mass of gas, which is conserved through it. ‘The amount of dissociation 
which ts going to take place behind the shock then depends solely on the 
physical constants of the gas itself and this energy. In particular, only the 
ratio of this energy to the energy required to dissociate the gas and the ratio 


of a characteristic density of dissociation pp to the free stream density py 
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appear. ‘lhese constants will be referred to as the ‘ equilibrium constants’ 
since their values can be determined from equilibrium considerations alone. 
The rate at which equilibrium behind the shock wave is achieved is 
determined by a constant C. The value (p, C)~! gives the time scale of 
the dissociation process behind the shock wave. 

In a similar way, the changes normal to the bow shock wave of a bluff 
body take place according to the Rankine—Hugoniot relations. Momentum 
tangential to the shock wave, and hence tangential velocity, is conserved. 
Since the density behind the shock wave increases from a value 7p, to 
higher values due to dissociation, the approximation of large density behind 
the shock used by the author (1956) and others (Chester 1956; 
Ivey, Klunker & Bowen 1948) and called the ‘ Newtonian approximation ’ 
is especially relevant. By using such an approximation it is possible to 
predict some of the geometrical properties of the flow, such as bow shock 
wave and streamline positions. In theory it is also possible to predict the 
pressure variation corrected from the crude * Newtonian’ value. In practice 
such a derivation proves rather laborious. But as the present investigation 
is based on qualitative rather than quantitative assumptions about the gas 
itself, the theory is regarded as sufficiently justified if it can predict a 
qualitative picture of the flow. 

The relevant non-equilibrium parameter in this case is denoted by A 
and is the ratio of the body size to the length scale of the dissociation 
phenomenon. For the normal shock wave the time scale of the dissociation 
process is (Cp,)~! and hence the length scale of the whole flow is U/Cp, 
where U is the free-stream velocity. ‘Thus A = /Cp,/U where / is the 
characteristic length of the body. For small A the body is small compared 
with the distance dissociation requires and thus for A = 0 we have the 
undissociated flow of a perfect gas with a constant ratio of specific heats 
y= 4. For large A, however, the main dissociation takes place at the 


/ 


bow shock wave, and dissociative equilibrium is maintained throughout 
the region between the shock wave and the body. Quite generally, we 
can make further predictions about the flow between shock and body by 
considering a ‘local’ value of A which we will denote by A. The local 
length scale of the dissociation is u/Cpy, where u is the local velocity and 
hence A=/Cp,/u. Near the stagnation point, A is large and hence 
dissociative equilibrium is achieved in this region. Also, on the basis of 
the approximation of the present theory, changes in the velocity along 
the streamlines are negligible and hence, near the body, where the stream- 
lines all originate in the region near the stagnation point, the velocity is 
small. ‘Thus, near the body, A is also large and we have dissociative 
equilibrium. Proceeding outwards along the normal to the body surface 
therefore, we have dissociative equilibrium at the surface whence the 
amount of dissociation decreases until it is zero on the shock wave itself. 
From the theory, dissociation profiles between the shock and body can be 
computed, showing in detail the amount of dissociation across the ‘shock 


layer’ at various stations on the body. 
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‘The amount of dissociation between the shock and the body, since it 
increases the density of the fluid in this region, will tend to decrease the 
distance between shock and body due to the contraction of stream tubes. 
Schwarz & Eckermann (1956) infer the amount of excitation of the 
vibrational modes of the molecules of polyatomic gases from measurements 
of the ‘stand-off’ distance of the shock wave from the stagnation point. 
This adds interest to the discussion (§4) of the variation of this length 
with the rate of dissociation. 


2. RATES OF DISSOCIATION 
The ideal dissociating gas has been discussed by Lighthill (1957) and 
the laws governing its equilibrium behaviour deduced. When thermo- 
dynamic equilibrium is not attained, however, it is necessary to consider 
the actual processes of molecular dissociation and atomic recombination. 
Equilibrium is achieved when the rate of dissociation of the molecules 
into atoms is equal to the rate of production of new moiecules by the 
recombination of atoms. ‘Thus, the net rate of dissociation 
dz. 
> ryn—lr 
where ~ is the ratio by weight of atoms dissociated to weight of atoms and 
and 7, denote the rates of dissociation and 
recombination respectively. In equilibrium, r, =r, and hence if we 
know either r, or rp, using the equilibrium theory (Lighthill 1957), we 


(2.1) 


molecules at a point and r, 


) 


) 
can deduce a relation between r,, and r, for the ideal dissociating gas. In 


particular we shall attempt to find a value for r,. Now dissociation of a gas 
molecule is accomplished by increasing the energy in the internal degrees 
of freedom to a point which 1s sufficient to overcome the binding forces 
holding the atoms together. ‘The necessary increase in internal energy is 
achieved by transfer from other forms of energy during a collision. If the 
collisions are sufficiently violent dissociation ensues. Hence, the rate of 
dissociation 7, will be to a first approximation proportional to the number 
of binary collisions of a molecule with another molecule or free atom such 
that the total energy available in the collision is sufficient to cause 
dissociation. We can therefore write 


rp = Cy(x, T)p(1—a)e- 2 *?, (2.2) 


where p is the density, 7 the temperature, D is the energy of dissociation 
and & is Boltzmann’s constant. C,(z, 7’) is a function yet to be determined. 
The factor p(1—«) gives the number of molecules per unit volume of the 
gas. Since we are primarily interested in binary collisions C, will, in general, 
consist of a sum of two terms, one corresponding to collisions between the 
molecules themselves and another for the collisions between atoms and 
molecules. At constant temperature the former will be proportional to 
(1—«) and the latter to x. ‘The factors of proportionality will depend on 
the respective collision cross-sections. An attempt to evaluate these has 
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been made by Wood (1956) and also Evans (1956) on the basis of simple 
kinetic theory. In view of the uncertainty attached to these values however 
the assumption will be made in this paper that C, is independent of x. The 
main variation in C,, it is thought, will come from the dependence on 
temperature, which we shall assume to follow an inverse-power law. Such 
a Variation is expected since the energy in many degrees of freedom (the 
rotational and vibrational energy of each colliding molecule, and their 
relative translational energy) may combine to make dissociation possible 
if their total exceeds D. ‘The inference of an inverse-power law in such a 
case (Hinshelwood 1940) rests on the fact that the proportion of states of 
a system in which the energy, made up of » ‘independent square terms’, 
exceeds a value D, large compared with 3nkT, 1s 
e-DkT(D/kT)*"—-! (An—-1)!. 
However, the precise value of 7 which ts appropriate is unknown, although 
it will be shown not to be crucial. The value » = 7, used below, includes 
the relative translational energy but the rotational and vibrational energies 
of only one colliding molecule, on the grounds that the vibrational energies 
are not fully excited, that the molecules sometimes collide with free atoms, 
and that it is uncertain whether any mode of collision could use all the 
available energy. We therefore assume the expression for the rate of 
dissociation to be 
r, = C,(1—a)T~* Dk? (2.3) 
where C and s are constants. ‘The dependence of the results obtained 
from this formula on the exponent s can be exhibited by choosing different 
values for s in the computation (in practice, 0 and 2-5). 
From equilibrium theory (Lighthill 1956), we note that for the ideal 


dissociating gas 





a” = PD emt (2.4) 


in conditions of thermodynamic equilibrium, where pp, is a constant 

characteristic of the gas. This is therefore the condition r, = rp. Hence 

an expression for rp of the form 

i, = ei pra? (2. 
PD 

will satisfy equilibrium conditions. Further, an expression of this form 

might be expected since recombination requires a three-body collision of 


wn 





) 


either three atoms or two atoms and a molecule and occurs, therefore, at 
a rate proportional to the square of the density px of free atoms. ‘The net 


rate of dissociation therefore becomes 


dy _ CTI (1—a)e-Det — & g2\. (2.6) 


dt \ Pp 


Comparison of expression (2.6) with the results of Wood (1956) and Evans 
(1956) shows that the main features agree. In these two papers different 


a 
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expressions of the dissociation rates are deduced. Ditterent values of s are 
obtained and also the variation of C with z is included. ‘The latter variation 
is however simply linear. Recombination is neglected in most of the work 
of Wood and Evans. 

The relation (2.6) will be used throughout this paper. ‘he variations 
of C from its constant value are probably small when compared with the 
variation of the exponential factor in (2.6) for any cases of interest. ‘This, it 
is hoped to predict a qualitative picture of the flow in non-equilibrium 


thermodynamic conditions. 


3. FLOW THROUGH A NORMAL SHOCK WAVE 


(1) Equations of one-dimensional motion 

In the region behind the Rankine—Hugoniot shock wave, the rate of 
dissociation is given by (2.6) with d/dt as the total rate of change with time. 
Since the flow is steady, however, the changes are purely convective and 
we can write (2.6) in the form 

— Cet ‘L(1—x)eP eT —g2l | (3.1) 
dx Pp) 

where x is the coordinate measured normal to the shock wave in a down- 
stream direction and wu is the velocity of the gas in this direction. 

In the region behind the Rankine—Hugoniot shock wave the gas is 
assumed inviscid and without heat conduction. ‘The equation of state for 
the ideal dissociating gas is then 


P— RT +0). (3.2) 
p 
Continuity of mass and momentum for strong shock waves requires 
pu = p, U, | 
p+ pu = p, U?,| 
where p is the pressure. ‘The constant on the right-hand side of these 
equations takes the value in front of the Rankine-Hugoniot shock since 
these quantities are conserved through it too. Finally the total energy 


— 
ww 
w 

— 


of the gas is conserved and thus 

t+ hu? = 1U?, (3.4) 
where 7 is the enthalpy of the gas. The five equations (3.1)—(3.4) then 
determine the problem completely provided that 7 is known in terms of 
the other variables. For an ideal dissociating gas the internal energy e per 
unit mass is of the form 


where m is the atomic mass. ‘The first term is the contribution to the 
energy from the various degrees of freedom (R is the gas constant for the 


undissociated gas) which is the same for both atoms and molecules. ‘This 
arises from the fact that although the atoms have only three degrees of 
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freedom, their mass is only half that of the molecules. Also the molecules 
of an ideal dissociating gas have only halt the energy in their vibrational 
degrees of freedom (and hence have effectively six degrees of freedom). 
‘The second term is the energy absorbed by dissociation. By use of (3.2) 
and (3.5), equation (3.4) yields 


(440)RT + lu2+ ~ = 10. (3.6) 


Furthermore, equations (3.2), (3.3) and (3.6) give 


5-2 (te)- J(- Gemeente, an 


OTe Teel ere PO 


where a new variable ». = mU?/D (the ratio of the kinetic energy of 
undisturbed gas to its dissociation energy) has been introduced. Equations 
(3.7) and (3.8) together with (3.1) then give 


dx : ‘ 
, F(x) (3.9) 

















for a certain complicated function F(«), and the variation of « with x is 


given by 





=f a. (3.10) 
Jo F(a) 
The thickness of the dissociation region can conveniently be defined as 
0-950, dy 
J0 F(a) ( ) 


where «, is the equilibrium value of the dissociation obtained from F(«) = 0, 
this being the distance required to reach 95°,, of the equilibrium value. 
Once the function x = x(x) is known from (3.10), it is then possible to 
obtain, by direct substitution, the values of all the other physical variables 
from (3.7), (3.8) and (3.2). 

Before the integral can be evaluated, however, it is necessary to determine 
the range of values of interest for the various parameters. The equilibrium 
constants D and p, are known and hence typical values of and pp/pp can 
be chosen. ‘The values» = 1 and 3, and pp/p,) = 10 and 10”, are used in the 
computations. The latter values are typical of the part of the atmosphere 
where the density is from a tenth to a hundredth of its sea-level value. Due 
to the exponential factor in (2.6) it is effectively only the logarithm of pp/pp 
which is important. 

The values of C and s are concerned with essentially the non-equilibrium 
conditions. C, which is unknown, determines the length scale of the 
dissociation process. Conversely, the thickness of the dissociation region, 
as observed in any experimental work, would determine the value of C 
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which could then be used in further computations—as those in the latter 
part of this paper. ‘I'he parameter s is given the values 0 and 2-5 as discussed 


above. 


(11) Computed results 

In the figures 1 to 4 various combinations of the parameters p, p,/po 
and s are used in the integration of (3.11) and the results are plotted 
for the ditferent flow variables against a non-dimensional length scale 
v, = (xp) C/R*°U'~*s). Table 1 is designed to give a key to the curves 





Curve pA PolPD s e A 
ees ® mS | 
(a) 1 10-* | 0 | 0-591 | 51 | 
| (b) } 1o-* | O | 0-204 | 570 | 
(c) 1 10-7 | O | 0-645 | 106 
(d) 1 w-- | O | 0-591 | => Wy 
(e) | 1 10-6 | 25 | 0-591 0-02 
(f) 1 io | 28 0-591 | 0-85 





Table 1. Key to figures 1 to 6. 


in figures 1 to 6. In figures 1 and 2, s = 0 and the variation of x (figure 1) 
and the density and temperature (figure 2) are exhibited for (a) u = 1, 
Po Pp LO=*s (b) LL 5, Po Pp 10%: (c) be I Po Pp 10>". In 
figure 3, s = 2:5 and the other variables are as in (a) and (4). The most 
characteristic feature of all these curves is the rapid variation near the 
Rankine—Hugoniot shock wave and the much slower variation near 
equilibrium. ‘This is almost entirely due to the exponential form of the 
energy distribution in the gas, as can be shown by replacing all the other 
factors by suitable constant values and repeating the integration. ‘This 
is most conveniently done by considering the temperature variation. ‘This 
procedure was adopted for result (a) and in figure 2 (4) the result obtained 
is plotted as curve (d). Curve (d) represents 


—+ = —0-0952 {e-1™ — e-l Tre) (3.12) 


dX, 
where 7 = 2R7T/U? and T,, is the value of 7, for equilibrium under the 
conditions of (a). 

An estimate of the distance required to reach equilibrium can be obtained 
by considering the distance required for the dissociation to reach 95°, of 
its equilibrium value. This value is denoted by A. A change in the total 
energy of the gas from the dissociation energy (u = 1) to half that value 
(4 = }) causes an increase in A by a factor of the order of 10. An interesting 
feature of the flow in this region is the only slight variation of the pressure 
there (figure 4). 
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‘The effect of variation of the parameter s can best be seen in figure 5, 
where the ratio of x to its equilibrium value ~, is plotted against v,/A. It 
will be seen that the increase in s tends to decrease the rapidity of the rise 
in x In order to contrast the difference of the results of the present 
investigation with the results observed in the absence of chemical changes 
the local length scale (x—x,)/(dz/dx,) is plotted in figure 6. ‘The rapid 
variation just behind the shock wave is again evident, whereas such length 
scales are approximately constant in vibrational relaxation behind shock 
waves. 
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Figure 1. Amount of dissociation x behind a normal shock wave (s = 0). 
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Figure 2 (a). The variation of density behind a normal shock wave (s = 0). 
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Figure 2(6). The variation of temperature behind a normal shock wave (s 
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Figure 3. Amount of dissociation 1 behind a normal shock wave (s = 2:5). 
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Figure 4. Pressure variation behind a normal shock wave (s = 0). 
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Figure 5. Comparison of the shapes of the dissociation curves. 
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Figure 6. The length scale of the dissociation process behind a normal shock wave 


4. FLOW PAST BLUFF BODIES 

Having discussed the behaviour of flow through a normal shock wave 
in a uniform flow, we can now proceed further and consider the flow behind 
the bow shock wave of a bluff body. In the previous work it was noted 
that the density behind a strong shock wave was large and that the effect of 
dissociation was to increase it still further. he theory developed by 
Busemann (1933, pp. 275-277), Ivey et al. (1948), and re-derived by the 
author (1956) makes use of the approximation of y near unity for 
a perfect gas with constant specific heats. Alternatively, this may be stated 
as requiring the density to be large—and, in general, this is the more correct 
statement. Since, therefore, in the dissociating gas this condition is satisfied 
more nearly than for a perfect gas, this type of approximation would seem 
to be a natural one to adopt. If therefore we use boundary layer coordinates 
x and y along and perpendicular to the body, the equations of motion, with 
viscosity and heat conduction neglected, are 


udu Ou uv. 1 op 0 
hox Ov or. phex '| 
| 
udv dv _ uw ldap 0. | 
hex Oy or. poey 


=~ 


0 0 
— (puk)+ — (pvhk) = 0, 
: i 

Di 1 Dp 

Dt p Dt j 
where u and v are the components of velocity in the x- and y-directions 
respectively, # dv and k dz are the elements of length in the v- and s-directions 
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respectively and D/Dt = (u/h)(é/ex)+7(é/ey). For an ideal dissociating 
gas (Lighthill 1957), 
; ae. 
i= (4¢+a)RT+ —x. (4.2) 
2m 
With use of the approximation (as in Freeman 1956) that the density is 
large throughout the region, (4.1) becomes 


enn eer a (4.3) 
v \ 
1a u> 
ge ae oe 4.4 
poy r ) 
s (puk) + 3 (pvhk) = 0, (4.5) 
Ox oy 
01 Cr 
—— Oo = 4.6 
, Ox ov sia, 


since changes between the shock and the body are large compared with 
those along the body. In terms of the stream function %, for which 
puk = cys/oy, puhk = — dys/ex, we have 


Op iu 


ay oe i= I(). (4.7) 


u = u(y), 
If, as assumed in the previous paper, the streamlines and shock lie close 
to the body, (4.7) gives p (the Newtonian plus centrifugal value as in 
Freeman 1956) in terms x and %. From (4.2), 
5 


on | 1 
RT = {1(u) - = al =: (4.8) 


and, by the equation of state of the gas (3.2), 
0 = p(x, ys) (; + O 4 . 
Tb) = (D/2m)x} \T $a)’ (4.9) 








where /(s) is the value of the enthalpy on the streamline % = const. 
Finally, using the equation (2.6) for the rate of dissociation, we obtain 





o% _ CphT~* fa —a)e- Dat P gah (4.10) 
<u) Pad 


Ww 


where h = h(x) in this approximation. The equation (4.10) together with 
(4.7), (4.8) and (4.9) becomes a first-order differential equation which 
can be solved along each streamline, i.e. for each value of %. The values 
of u(y) and J() can be deduced from their values on the shock wave where 
the streamline ¢% crosses it. For this purpose, the shock wave can be assumed 
to lie along the body itself. 
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Having in this way obtained the amount of dissociation along each 
streamline we then find it possible to obtain the positions of the shock wave 
and streamlines from the formula 


y= | dp (4.11) 


for x constant. ‘This may be written 
% 47(b)—(D/2m)x' (1 +) 
| » R(x). p(x, b). (4+a).u(eb) 


which can be evaluated when the solutions of (4.10) are known. ‘The theory 
will again break down in regions where p becomes small, or p becomes 





dis, (4.12) 


small, as noted in the previous paper (Freeman 1956). This occurs at 
some point on most bodies, since the centrifugal forces tend to decrease 
the pressure at the body surface. 

Near the body, the velocity u is small due to the fact that the streamlines 
there originate at the stagnation point. In the case of two-dimensional 
tow this requires a higher approximation to u to be found to make (4.12) 
converge, although this is not necessary (see Freeman 1956) in the three- 
dimensional case. Since the condition u — 0 at the body requires that 
the terms in the curly bracket of (4.10) should vanish, we infer that 
dissociative equilibrium is achieved in this region. ‘Thus along the body 
itself the amount of dissociation is determined by 


2 »>—-D kT 
a (4.13) 


l—« p 





Physically this means that the time scale of the flow is so large that the gas 
has time to adjust itself to equilibrium. If this is the case, it is possible to 
obtain x along ys = 0, as the variation of p(x,) is already known. Hence 
it is relatively easy to obtain a second approximation to uw near the wall by 
retaining the pressure term in the first equation of (4.1) to a first 
approximation and putting ys = 0. Hence 


u = u(ys)— [" om antl, (4.14) 


where p is determined by (4.9). It is possible to develop the theory for the 
two-dimensional flow provided the above higher approximation is 
introduced. 

In the following, however, we concern ourselves solely with axially 
symmetric blutf bodies since these are likely to be the most important 
from a practical viewpoint. Further, we shall consider the particular case 
of a sphere. ‘This is the simplest type of body and will give the behaviour 
in regions near the stagnation point correctly for all bluff bodies with the 
same local radius of curvature. The approximation does however break 
down (Freeman 1956) before the streamlines reach an angle of 60° from 
the front stagnation point due to the decrease in pressure along the body. 
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For a sphere of radius a, we introduce the new variable € defined by 
y = hp) Ua?sin?é. This gives u = Using, I = }U*cos*é and 


fsin 34+ sin® g| " . 
= { ———— >» U? 4.15) 
P | 3sin@ | ™ ( 

(see Freeman 1956), where @ is the angle measured from the front stagnation 
point. Equation (4.10) may then be written in the form 


Om 4+a J 
20 (+a) —[xia))| 
_ 2pof sin 36 + sin? (4+ %)x? 

pp \ 3sin@cos*é /(1+«)(1— [a/p;])) 


2aCpo(" {_sin3é+sin?é | 1, 
U }. \3sindsinécostéf 


Here, we have assumed that s = 0 and written p, = wcos*é. The latter 
integral can be evaluated to give 





(1—a)e“@+91-a) — 


(4.16) 





where 


® = 





2aCp, [6 —€ + sin 26 — sin 2é + sin’ € log(tan $0/tan $€)] 

U 3 sin € cos? € 
Substituting in (4.16) the value of @ from (4.17), we obtain a differential 
equation in x and © along the streamline € = constant. This must be 
solved numerically however since (4.17) gives only an implicit relation 
for @ and equation (4.16) is extremely complex. When (4.16) is solved, 
the position of the streamlines may be obtained from 
cos (1 + )(1— fala 
Jy (sin30+sin?é)(4+«) ’ 
and the position of the shock is then derived by putting € = @ in (4.18). 
Near the front stagnation point, it is possible to extend the theory a little 
further. In this region, since sin€ ~é and cosé~1, (4.16) may be written 


QO = (4.17) 





(4.18) 





bol Ww 


bs = 




















0x 4+ f 2p a%(4+a)  ) 
= = 4 (1—a)e—4@+o—-o) _ 0 \, (4.19) 
20 ~ (Tray —[alepi' pp (+ay(1—[xn))f 
where ve (=) = 2aC po 
é l 

and (4.18) simply becomes 

y_ 1 pf (+a) = fain) 

-= = 1E. (4.20 

a 20 J, (4+) “~ ) 
Assuming that « is a function of © only, we then have 
dz (4+) ( 2p 4°(4 +) 
— = 4(1—a)e—4+ve-—o __ = », (4.21) 
do) (isa\i-fepi pp (1a) [xe 
and from (4.20) 

~  - [" (1+«)(1—[a/n]) do (4.22) 


a ; ZN . A(@—&) € (4+)(1+[0/A])? , 
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‘The stand-off distance, 5, is then given by putting 6 = € in (4.22), and hence 
we obtain 





oY l [ (1+)(1—[z/u]) do 
2A J, (4+2)(1+[0/A])?? 


where 


2p, a2(4+% oa 
x)e-4+a)u—a) — Po ( ) —* 


pp (1+2)(1— [z/]) 








5. NUMERICAL COMPUTATIONS 

In the preceding section, a theory has been developed to give the flow 
pattern for a uniform stream of an ideal dissociating gas past a sphere. 
In view of the complex nature of the equations, however, it was not possible 
to state explicitly any results. It is necessary therefore to continue by 
adopting numerical methods and restricting consideration to particular 
cases. In view of the fact that probably most of the variation is determined 
by the exponential factor, it was decided to take one particular flow (with 
= 1, po/pp = 10-*) and attempt to exhibit the dependence on the 
coefficient A which determines the effective time scale of the dissociative 
process. The results obtained are exhibited in figures 7 to 11. Figure 7 
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Figure 7. The amount of dissociation along the streamline € = 5° for various values 
of A in flow past a sphere (yu 1, po/pp = 10-*). 
. shows a typical result for the dissociation along a streamline (€ = 5 ). 


It will be noticed that the variation along each streamline is similar to that 
behind a normal shock wave. The rate at which dissociation takes place 
however depends very much on the value of the time scale factor A. For 
A = 0, we obtain the solution for a perfect gas with y = 4, whereas for 
large we approach the solution for a gas in dissociative equilibrium. 
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Figure 8. The dissociation profiles between the shock wave and the surface of the 
sphere (ue 1, Po Pp = 10-8), 
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Figure 9. The amount of dissociation « along the dividing streamline between the 
shock wave and the stagnation point for various values of A(u = 1, po/pp = 107°). 
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Figure 10. The shock shape for various values of A(u = 1, po/pp = 10~°). 
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Figure 11. The ‘ stand-off’ distance 6 as a function of A(u = 1, po/pp = 10~*). 


Figure 8 shows the dissociation profiles at @ = 10° and 20°. he dissociation 
begins to be important closer and closer to the shock wave (€ = @) as A is 
increased. For the smaller values of A (< 1) the dissociation takes place 
almost completely in the region near the body where the velocities are small. 
In figure 9, the results obtained at the stagnation point from equations (4.21) 
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and (4.22) are plotted to show the variation along the dividing streamline 
ot x, the amount of dissociation. The rapid rise of « in the regions of low 
velocity when A is small is again evident. From a physical viewpoint 
however, figures 10 and 11 probably give the most useful information. 
The distance of the shock wave from the body can be easily measured from 
Schlieren or shadow photographs and consequently any variation with 
dissociation observed. In figure 10, the distance of the shock wave away 
from the body for various values of A is shown, and in figure 11 the stand-off 
distance itself is plotted as a function of A. It will be seen that between 
the limits of a gas with constant specific heats and a dissociated gas in 
thermodynamic equilibrium the distance of the shock wave away from the 
body can change by as much as 50°. It is hoped that these results, although 
only of limited application since they are obtained in this one particular 
case, are of value in showing that appreciable changes in the geometry of the 
flow field may be expected. 


The author wishes to express his sincere thanks to Professor M. J. 
Lighthill, F.R.S., for his help and encouragement during the preparation 
of this paper. This paper is published by permission of The Director, 
The National Physical Laboratory. 
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SUMMARY 

Consideration of the structure of wind-generated waves when 
the duration and fetch of the wind are large suggests that the 
smaller-scale components of the wave field may be in a condition 
of statistical equilibrium determined by the requirements for 
attachment of the crests of the waves. A dimensional analysis, 
based upon the idea of an equilibrium range in the wave spectrum, 
shows that for large values of the frequency w, the spectrum ®(w) 
is of the form 

D(w) ~ ag" > 

where x is an absolute constant. ‘lhe instantaneous spatial 
spectrum ‘i‘(k) is proportional to k~* for large wave numbers , 
which is consistent with the observed occurrence of sharp crests 
in a well-developed sea, and the loss of energy from the wave 
system to turbulence and heat is proportional to p,,u,, where p,, 
is the water density and u, the friction velocity of the wind at 
the surface. This prediction of the form of @(w) for large w 
with x = 7-4 10-3, agrees well with measurements made by 


Burling (1955). 


1. ‘(THE CONCEPT OF AN EQUILIBRIUM RANGE 

When a turbulent wind commences to blow across a large sheet of 
water initially at rest, a wave system is initiated and develops under the 
continued action of the wind. According to a theory described in a recent 
paper (Phillips 1957), a considerable part in this process is played by a 
resonance between the convected pressure fluctuations on the surface 
arising from the turbulence in the wind and the modes of free surface wave 
propagation, resulting in a rapid growth of the wave amplitudes. This 
theory, being concerned with the early stages of wave growth, is based 
upon a linearized surface boundary condition which neglects terms whose 
magnitudes relative to those retained are of the order of the wave slope. 
It is predicted that, after a short initial stage, the root-mean-square amplitude 
of each Fourier component of the wave system increases with time ¢ as t!?. 
Now, if the turbulent wind continues to blow for a sufficiently long time 
and if the fetch is sufficiently large, the wave amplitudes continue to increase, 
so that ultimately the linearized approximation of this theory becomes 


inadequate and interactions among the Fourier components of the wave 
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field become increasingly important. This present note is concerned with 
situations in which the fetch and duration of the wind are such that the 
amplitude of at least some components of the wave field can no longer 
be considered as ‘infinitesimal’ and interactions among these components 
are appreciable. 

Some valuable qualitative information on the properties of the wave 
system can be inferred from even casual observations of wind-generated 
waves on the ocean. In the first place, the surface displacement at a given 
point is a random function of time, or at a given instant is a random function 
of position in the sense that the detailed configuration of the surface is not 
reproduced by apparently identical meteorological and other conditions. 
This situation is similar to that found in the study of turbulent motion and 
suggests that similar statistical concepts should be applied. In particular, 
only the average properties of the waves can be regarded as significant 
experimentally or predictable by any reasonable theory. Secondly, it is 
commonly observed that in a well-developed sea, or one in which there 
are waves of ‘finite height’, a characteristic property is the occurrence of 
fairly sharp wave crests with intermittent patches of foaming (‘ white-caps ’ 
or ‘white horses’) which seem to develop when the crests can no longer 
maintain their attachment to the remainder of the water. 

The exact nature of this detachment has been considered by many 
people, and perhaps the most pertinent investigation is that of G. I. ‘Taylor 
(1953) who studied the limiting configuration of standing waves in a wave 
tank. In the open ocean, in deep water, it is likely that the limiting 
configuration of part of the water surface is reached when the local downward 
acceleration of the fluid near the wave crest is equal to g, the acceleration 
due to gravity. Ifthe surface configuration changes in such a way that the 
local acceleration tends to increase beyond this value, that is, in such a way 
as to put the fluid near the crest in a state of tension, (as happens if the wave 
crests tend to become sharper), the surface breaks and the crest of the wave 
detaches. The formation of these sharp crests leading to detachment is 
observed when two waves of the physical wave pattern run together, or 
when a wave, travelling with its phase velocity, moves into a region where 
the energy density of the wave field is locally high. A typical sequence 
of events in time can be summarized by saying that a wave which is initially 
rounded, develops a sharp crest which may either subside or detach and 
disappear, dissolving amidst a foaming patch of turbulence and leaving 
the wave again rounded. An alternative way of visualizing the same 
process is to consider an instantaneous photograph of a large region of the 
sea surface. Such a photograph would show a large number of rounded 
wave maxima where the local surface acceleration was less than the 
gravitational acceleration g, and also a smaller number of sharp wave 
crests where the surface acceleration happened to be exactly equal to g, so 
that at these crests the surface was on the point of breaking. It would 
also show some white foam patches marking the demise of sharp wave crests 
that had reached the limiting state at previous instants and broken. 
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We now enquire whether these observations can be translated into 
mathematical language concerning the instantaneous spectrum of the 
surface displacement. ‘The occurrence of scattered sharp wave crests as 
a transient limiting configuration corresponds to the occurrence of dis- 
continuities of surface gradient which in spectral terms corresponds to the 
existence of a certain form of the spectrum of high wave-number. ‘The 
properties of the instantaneous spatial spectrum at these high wave-numbers 
will therefore be determined by the physical parameters that determine 
the extreme configuration of the surface in the limiting condition, the 
particular property that is relevant being the magnitude of the discontinuity 
in surface slope developed. It seems likely, therefore, that in a well- 
developed sea there exists a range of large wave-numbers over which the 
wave spectrum is statistically determined by the physical quantities 
governing the conditions for attachment of the wave crests. Similar 
remarks can be made concerning the frequency spectrum of the surface 
displacement at a given point, where rapid changes in the surface displace- 
ment are associated with the movement past the point of observation of 
occasional sharp wave crests near the limiting configuration. 

‘The basic hypothesis of this paper can therefore be stated briefly as 
follows. In a well-developed sea, generated by the wind, there is an 
‘equilibrium range’ of large wave-numbers (or high frequencies) in the 
spectrum, determined by the physical parameters that govern the continuity 
of the wave surface. 

There are two remarks that should be made at this point. ‘The first 
is that the concept of an equilibrium range is suggested by consideration 
of the asymptotic form of the spectrum for large wave-numbers (or 
frequencies) and that we have little direct evidence upon which to base 
an a priori estimate of the smallest wave-number to which we would expect 
it to be applicable. It is clear that a necessary condition for the existence 
of an equilibrium range over a certain part of the spectrum is the existence 
of appreciable non-linear interactions among these wave-numbers, but it 
is not clear whether this condition alone is sufficient. The results of some 
measurements described in $3 offer good a posteriori evidence that it may 
indeed be sufficient but it may be difficult to justify such an assertion in 
advance. ‘The second remark is that the magnitude of the wave spectrum 
in the equilibrium range represents an upper limit, dictated by the require- 
ments of crest attachment. In the early stages of wave generation, the 
equilibrium value may not have been attained at any point in the spectrum, 
although the wave slopes may be such that non-linear interactions are not 
negligible, and in a decaying wave system the wave spectrum over the 
relevant range may have fallen from its equilibrium value through the 
damping of the components of shorter wavelength. Fruitful consideration 
of these situations is difficult, even on simple dimensional grounds, because 
of the additional parameters involved, and we will restrict our attention to 
the equilibrium range, which might be expected to be attained when a 
statistically steady wind of sufficient duration and fetch continues to act 


upon the water surface and supply energy to the wave motion. 
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2. DIMENSIONAL ANALYSIS 

Since we have defined the equilibrium range of the wave spectrum to 
be that part governed by the requirements of attachment to the wave crests, 
it would be expected to be independent of the fetch and duration of the 
generating wind. ‘The physical quantities to be considered in a dimensional 
analysis, then, are the densities p,, and p,, of the air and the water, the friction 
velocity u, representative of the wind speed, the surface roughness length 
2», the gravitational acceleration g, the surface tension and viscosity of the 
water 7 and v, together with the wave-number & or the frequency w. But 
not all of these parameters will be important. Firstly, provided we restrict 
attention to wave-numbers and frequencies well below those associated 
with capillary ripples, so that 

; ) g3) 1/4 


) 1/2 
Pus < 4p wi" (1) 


or 
ao Ce ea 

the influence of surface tension is unimportant. Furthermore, the direct 
effect of viscosity is manifest in wave damping, the logarithmic decrement 
being y = 2vk? (Lamb 1932, §§ 348, 349), and, under the conditions with 
which we are concerned, the rate of energy transfer from the wind to the 
waves is such that this damping can reasonably be ignored. The condition 
(1), in the present situation, suffices to ensure that the influence of viscosity 
in the water is also unimportant. ‘Thirdly, as soon as we agree that T can be 
omitted from consideration, it will be observed that the dimension mass 
occure only in p, and p,,, so that if we are to be concerned with quantities 
that do not contain this dimension, then the densities can occur only as 
functions of the combination p,,/p,,.. Under oceanographic conditions, with 
an air—water interface, this ratio is virtually constant, although under more 
general conditions, at the interface between two fluids, it may vary between 
zero and unity and its variation may demand consideration. 

Of the remaining parameters mentioned above, there is some evidence 
that the surface roughness length z, is not independent of the others. One 
might reasonably expect the surface roughness to depend upon the wind 
speed, specified by u,, and Ellison (1956) suggests that 

Bp K us /g. (2) 
Accurate measurements of z, over the sea are difficult to make and the results 
of different workers are by no means consistent. However, Ellison quotes 
the results of some careful measurements by Hay (1956) which support 
a relation of this type and further evidence has been collected by Charnock 
(1955). The length z, is determined by the height of the short steep waves 
on the water surface and is probably of the order of one-thirtieth of the 
height of these waves (judging from experience with solid roughness 
elements), so that Hay’s measurements, made for values of z, less than 
():3 cm, probably to refer to conditions under which the amplitude of such 
waves is less than about 10cm. However, if the expression (2) is valid 
for these relatively short finite amplitude waves, it is likely to be true for 
larger values of z, and u,, and we can reasonably conclude that the equilibrium 


kR<€ 
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range properties of the wave field are determined only by the parameters u, 
and g (and in the general case p,/p,,.), together with k or w, provided the 
wave-numbers and frequencies satisfy the conditions (1). 

‘he functional form of the wave spectrum in the equilibrium range can 
now be found readily. ‘The frequency spectrum of the surface displacement 
&(x,?) at a fixed point is defined as 


or 








@(w) = (27) | &(x, t)€(x, t+ 7)e—" dr, (3) 
and the instantaneous wave-number spectrum as 
V(k) = (27) | &(x, t)€(x +4, the" dr, (4) 


the integration being over the whole surface. ‘The function ®(w) has 
dimensions L?7 and '’(k) has dimensions L*. As we have seen, the form 
of these expressions for large values of their arguments is determined by 
the limiting configurations attained by the surface under the requirements 
of attachment of the wave crests. The geometry of the limiting shape near 
the sharp crests is determined* by the condition that the downward 
acceleration should not exceed g, so that the asymptotic forms of (3) and (4) 
would not be expected to involve u,. An increase in u, would have the 
effect of increasing the rate at which wave crests are passing through the 
transient limiting condition of an abrupt change in surface gradient but 
should not influence the geometry of such a sharp crest itself. 

As far as the functions (3) and (4) are concerned, then, the relevant 
physical parameter seems to be simply g, as well as either w or k as the case 
may be, so that it follows immediately on dimensional grounds that 


D(w) ~ «gw, (5) 
where ~ is a constant, and 
V(k) ~ f(@)k-4, (6) 


where @ is an angle specifying the direction of the vector wave-number k 
and f(#) is determined by the directional distribution of the smaller-scale 
components. These expressions might be expected to refer to ranges of 
frequency and wave-number such that 

wy Kw < (4p, °T 1) 4; Ry RK (pwgT)*, (7) 
where w, and ky represent the smallest frequency and wave-number of the 
wave field for which non-linear effects are important. 

It is interesting to notice that the expression (6) is consistent with the 
particular form that the limiting configuration assumes, that is, with the 
occurrence of sharp wave crests. It can be shown that if a function has 
simple discontinuities in gradient at one or more points, then its Fourier 
transform is asymptotically proportional to k~?. ‘Thus the wave spectrum, 
which is essentially the mean-square of the Fourier transform, is proportional 


* We here exclude a different possible type of surface instability in which the 
sharp crests may be ‘ blown off’ by very high winds. 
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to k-* for large k. Indeed, this remark suggests that (6) can be obtained 
alternatively by replacing our argument that uw, is not involved in the 
equilibrium range of the spectrum at any instant by the specific premise 
that the limiting configuration is characterized by the occurrence of sharp 
wave crests. Then, it follows from the statements immediately above that 
‘(k) « k-4, with a coefficient of proportionality that may now involve u, 
as well as g. But this coefficient must be dimensionless, and since there is 
no dimensionless combination of u, and g, the coefficient must be constant 
for each direction of the vector wave-number k. This argument is in some 
ways more attractive than the previous one leading to (6), though the specific 
nature of the premise used precludes an immediate deduction of the result (5) 
for D(w). 

Another important property of the wave field that is determined by the 
equilibrium range of the spectrum is €, the mean rate at which energy is 
lost from the wave motion on unit area of surface either to turbulence by 
wave breaking or by direct dissipation in the smaller scale components. 
The process of formation and detachment of sharp wave crests corresponds, 
though possibly not in a simple way, to energy transfer across the equilibrium 
range to high frequency components with its subsequent loss from the wave 
motion. If the range of frequencies from which energy is extracted from 
the wind is distinct from the range over which energy is lost, then « represents 
the rate at which energy is transferred across the equilibrium range of the 
wave spectrum, and so is determined by uw, and g alone. ‘The dimensions 
of « are (energy)/(time x area), so that 

«= Bo. ur (8) 
where § is a dimensionless constant. If the internal energy dissipation 
in the waves is neglected, « represents the difference between the rate of 
gain of energy from the wind and the rate of increase of energy of those 
components of the wave system that are still developing and have not yet 
attained a statistical equilibrium state. 


3. MEASUREMENTS OF THE HIGH FREQUENCY 
COMPONENTS OF THE WAVE SPECTRUM 

Among the most reliable measurements yet available of the wave 
spectrum (w) at large frequencies are those of Burling (1955), made on 
Staines Reservoir, Middlesex, England, using the capacitance wire recorder 
developed by ‘Tucker & Charnock (1955). This instrument records 
accurately the high frequency components of the wave system and so is 
particularly suited to an investigation of the existence and properties of 
any equilibrium range in the wave spectrum. A considerable number of 
wave records were made over fetches of from 400 m to 1350 m, with wind 
velocities at a height of 10 m ranging from 500 to 800 cm/sec. The results 
of taking Fourier analyses of these records are shown in figure 1. When the 
frequency w is small, the spectrum curves depend quite strongly on the 
fetch and meteorological conditions, but their most remarkable property 
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is that when w is large, the curves obtained under different conditions very 
nearly coincide and become apparently independent of the fetch and of 
the strength of the wind. So striking was this observation that Burling 
himself suggested the possibility of some type of equilibrium structure, 
but the idea was not followed up. He found empirically that the mean 
values of his observed spectra in this range obeyed a relation of the type 
P(w) = 7-0 x 108w* 
inc.g.s. units, where P(w) is as defined in (3) and the frequency w is expressed 
in radians per second. ‘This frequency dependence is the same as has been 
predicted by (5) and the observed value of the coefficient indicates that the 
dimensionless constant x in (5) is approximately 7-4 x 10-3. These measure- 
ments appear to offer good support both to the original hypothesis of the 
existence of an equilibrium range and to the result (5) based on this idea. 
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Figure 1. Spectra of wind-generated waves measured by Burling (1955). The 
cluster of lines on the left are representative of the spectra at low frequencies 
for which equilibrium has not been attained. On the right the curves merge 
over the equilibrium range, and the broken lines indicate the extreme measured 
values of @(w) at each frequency w. ‘The crosses represent the mean 
observed value at each w, and the heavy line the relation ®(w) = xg*w~° with 
ya 7 xD 


It is interesting to compare these results with those obtained by other 
oceanographers which are summarized in various empirical or semi-empirical 
spectra proposed as descriptions of the sea surface. Bretschneider (1958) 
proposes a spectrum of the form 

S(w) 427w 5 exp; —():675(27 Tw)*}, 
where 7 is a ‘mean wave period’ which depends upon wind speed and the 
state of development of the sea. When w is large, Bretschneider’s spectrum, 
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which is dimensionally sound, reduces to the form (5) derived above, and 
Bretschneider finds that the value x = 7-4 x 10-* obtained from Burling’s 
measurements is consistent with his own observations. On the other 
hand, Neumann (1954) proposed the expression 


I(w) = Cu exp} — 2g%/0® W%, 


for a ‘fully developed sea’, where W is the ‘wind speed’ and C is a 
dimensional constant equal to 1-40 x 10'cmsec~*, I(w) being normalized 
in the same way as our definition (3). The presence of this dimensional 
constant C cannot easily be justified on dimensional grounds, since it must 
be determined by some physical parameters yet is independent of w, W (or u,) 
and fetch, presumably leaving only g which cannot account for the 
dimensions of C. If we put this important objection aside, however, and 
examine Neumann’s spectrum as it stands, it is clear that when a is large 


I(w) ~ Cw, (9) 


which is independent of the wind speed and so contains one of the properties 
of the equilibrium range as described in this note. But the frequency 
dependence is as w~® rather than as w~*, which is required by dimensional 
analysis and supported by Burling’s observations; and the latter deserve 
more weight in view of the more adequate and precise instrumentation. 
Furthermore, the magnitude of /(w) predicted by (9) is only about 25°, 
of the values observed by Burling in this frequency range. However, 
in fairness to the Neumann spectrum as an empirical formula for describing 
ocean waves, it should be mentioned that the data from which Neumann 
obtained his spectrum cover frequencies much lower than those observed 
by Burling, so that (9) represents a considerable extrapolation from frequency 
ranges over which this spectrum has been found to give a reasonable empirical 
fit. Since it seems that the Neumann spectrum is unlikely to represent a 
genuine physical law on account of its dimensional inconsistency, its failure 
at these high frequencies is hardly surprising, and of course does not damage 
its usefulness in practical wave forecasting. 


This work was supported by Contract Nonr 248(38) with the Mechanics 
Branch of the Office of Naval Research. 
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REVIEWS 


Water Waves: The Mathematical Theory with Applications, by 
J. J. Stoker. New York: Interscience Publishers, 1957. 567 pp. 
$12.00 or 88s. 


Most students of fluid dynamics are familiar with the chapter on water 
waves in Lamb’s Hydrodynamics. Fewer, however, are aware of the many 
interesting developments in the theory of water waves since the last 
edition of Lamb’s treatise; developments stimulated no less by the 
temporary demands of the second World War and the constant need for 
engineering and oceanographic knowledge than by the intrinsic mathe- 
matical interest of the subject. Professor Stoker’s long-awaited book is 
the first serious attempt, since the war, to give a mathematical account of 
this rapidly growing field. 

Two things about the book stand out immediately. First it is written 
almost purely from the mathematical point cf view. ‘There 1s little of the 
quantitative comparison with observation so desirable in ali branches of 
fluid dynamics, and not least in the study of water waves. 

Secondly, only certain parts of the field are covered; for example, 
there is no account of surface tension or viscosity, of internal waves, or of 
standing waves of finite amplitude. At the same time, some very 
specialized topics, such as the wave pattern produced by a ship steering 
an exactly circular course, have been treated in detail. 

The book is divided into four parts. ‘The first part, consisting of two 
short chapters, is devoted to a derivation of the general equations of 
irrotational motion for a non-viscous fluid, and to some well-known 
theorems required later in the book. ‘The peculiar difficulty arising from 
the non-linear boundary conditions at the free surface is noted. ‘Two 
basic approximations are described, namely the * small-amplitude’ 
approximation of Airy and Stokes, where the solution is developed in powers 
of the ratio of wave height to wavelength; and the ‘ long-wave’ or 
‘ shallow-water’ approximation of Boussinesq and Rayleigh, which is more 
appropriate when the ratio of mean depth to wavelength becomes small. 

Part II of the book then describes some solutions with the small- 
amplitude theory (to a first approximation only) and Part III describes 
solutions with the shallow-water theory. Finally, Part IV contains two 
somewhat dissimilar problems, namely the problem of the breaking dam 
and the question of the convergence of the small-amplitude method of 
approximation. These are loosely associated as ‘* problems in which the 
free surface conditions are satisfied exactly ”’. 

In Part II a conventional treatment is first given of simple harmonic 
waves in water of uniform depth, with an interpretation of group-velocity 
both as the velocity of propagation of a wavelike disturbance (represented 
in the form of a sum or integral) and as the velocity of propagation of 
energy. Chapter 4 treats the motion due to a harmonic distribution of 
pressure over part of a free surface. One of the principal mathematical 
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tools used in this section of the book is complex variable theory, combined 
with suitable radiation conditions at infinity. In chapter 5 the small- 
amplitude theory of waves on sloping beaches is described using the 
artificial assumption of a logarithmic singularity at the shoreline. In 
reality, of course, the energy is dissipated in some other way, either in 
friction at the bottom or by breaking of the waves, so that non-linear 
terms would certainly be important. However, the linear problem has 
been elegantly solved by Peters and Roseau for arbitrary angles of inclina- 
tion of the beach. Professor Stoker gives also his own method of solution 
for special angles 7/2. Similar complex-variable methods are applied to 
the problem of waves breaking against a flat ‘ dock’ and against a vertical 
‘cliff’. In the latter case the solution is hardly realistic, since no account 
is taken of the reflected energy. The section includes a solution of 
Sommerfeld’s problem of the diffraction of waves round a vertical wedge. 

The second subdivision of Part II, entitled * Unsteady motions ”’, 
treats transient motions started from rest, using a Fourier transform 
technique combined with Kelvin’s. stationary phase approximation. 
Similar techniques are used in the third subdivision, dealing with waves on 
running streams (or waves caused by a moving pressure pulse). In the 
latter problem the question arises how to determine the uniqueness of 
the solution; for a free wave of arbitrary amplitude may always be added 
to any solution of the free surface conditions. One way is to assume a 
small but finite frictional dissipation of energy. Professor Stoker makes the 
more satisfactory assumption that the steady state is the result of starting 
the motion from rest and waiting until the transient motion has died down. 
According to this assumption, the transient energy is radiated away to 
infinity. 

‘The three-dimensional pattern caused by a moving pressure pulse 
(Kelvin’s ship-wave problem) is treated in chapter 8, in some detail, 
although with no discussion of the ‘ cusp-line’ where the wave form 1s 
expressible in terms of Airy functions. Part I] ends with a chapter giving 
a linearized theory of ship motion, in which it is also assumed that the ratio 
of the ship’s thickness to length is small. ‘Though similar in some respects 
to Haskind’s theory, it ignores terms that would lead to a damping of the 
pitching, surging and heaving motions, and so cannot account for these 
phenomena. As in all the preceeding chapters, no quantitative compari- 
son with observation is given. 

In Part II] we come to the author’s discussion of waves in shallow 
water. Chapter 10, which is largely a reproduction of a paper by the 
author (Comm. Pure Appl. Math. 1, 1948, 1), is the longest in the book and 
is the one whose validity seems most open to question. ‘The theoretical 
background is roughly as follows. 

In water of finite depth, the convergence of the small-amplitude method 
depends upon the parameter 


€ = akj(kh)s 


being small compared with unity, where a denotes the wave amplitude, 
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h the mean depth and & is the inverse of some typical horizontal dimension 
such as the wavelength. When the depth diminishes, so that € is of order 
unity or greater, another method of approximation (the shallow-water 
approximation) may be used; this involves expanding the solution in 
powers of (kh)°. 

It is easily shown that in the first approximation to the shallow-water 
theory all disturbances are propagated without change of form, with 
velocity (gh)'*. In the second approximation any progressive wave has 
in general a small second-order deformation; if we set this deformation 
equal to zero we obtain a non-linear equation for the first approximation, 
which yields a class of steady progressive waves, the * cnoidal’ waves of 
Korteweg & de Vries (1895). As the amplitude diminishes these tend to 
the well known sine-waves of Stokes, and as the wavelength is increased 
one obtains, in the limit, the solitary wave of Rayleigh and Boussinesq. 
Convergence both for the solitary wave and for cnoidal waves has recently 
been proved by Friedrichs & Hyers (1954) and by Littman (1957). ‘The 
work of Korteweg & de Vries showed further how a shallow-water wave 
of arbitrary form will gradually change shape, becoming sometimes 
steeper in front and sometimes less steep according to circumstances. 

Professor Stoker makes only a cursory reference to the work of 
Korteweg & de Vries, and takes no advantage of their excellent elucidation 
of the nature of long waves. His account of these matters, particularly 
of the limitations inherent in the shallow-water theory, is too diffuse and 
leaves much in doubt. In the method of treatment with which most of 
chapter 10 is concerned, all terms of order (kh)* are entirely neglected, 
giving rise to what may be called a zero-order theory. ‘The equations are 
indeed convenient to handle, being analogous to those governing the 
propagation of plane sound waves of finite amplitude; solutions can be 
constructed using the method of characteristic curves. However, this 
approach leads to the conclusion, for example, that any wave advancing into 
an undisturbed region must become steeper and eventually break, and that no 
steady wave can exist (apart from the trivial case of an unperturbed stream). 
‘This is directly contradicted by the existence of solitary and cnoidal waves 
which are well established both by theory and observation. 

The zero-order theory may of course have a limited region of applic- 
ability, and a careful discussion of this point is much to be desired. But, 
since no account is taken of vertical accelerations, the theory can hardly 
apply to breaking waves; Professor Stoker’s attempts to press the theory 
this far must therefore be viewed with reserve. His theoretical profiles 
of breaking waves are indeed unlike what is observed, as can be seen from 
Professor Munk’s photographs shown in the book. 

Chapter 11, entitled ‘‘ Mathematical Hydraulics’, is unlike the rest 
of the book in that mathematical elegance is sacrificed in an attempt to 
predict a concrete phenomenon, namely floods in the Ohio and Mississippi 
Rivers. An empirical law of turbulent friction is assumed, and good 
agreement is obtained with observation. ‘Uhis chapter follows closely two 
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reports to the U.S. Army Corps of Engineers by the author, in conjunction 
with Isaacson and ‘l’roesch. 

The final chapter of the book contains the two loosely associated 
problems mentioned earlier. ‘lhe first of these, the breaking of a dam, 
is solved by the author in Lagrangian coordinates, using a series expansion 
in powers of the time. In fact only the terms as far as & are evaluated. 
The problem appears to be similar to others (for example a collapsing 
mass of water with a sharp crest), in which difficulties are encountered in 
the terms in ¢4, and one would be happier to see the approximation carried 
out to terms of higher order. 

This chapter also contains a discussion of the existence of progressive 
waves of finite amplitude, and another version of Levi-Civita’s proof for 
the convergence of the small-amplitude approximation. I feel that undue 
emphasis is given to this existence proof while the significance of other work 


relating to the problem is neglected. In Levi-Civita’s work the ratio of 


wave-height to wavelength is limited to about 1/200 (see Hunt 1953) 
whereas experimentally much steeper waves are known to exist. ‘The 
surface profile of the * highest wave ’ (that is to say, a wave having a sharp 
crest with 120 angle) was calculated long ago by Michell (1593). Further, 


Havelock (1918) extended Michell’s method so as to compute the form of 


a wave of any amplitude less than the highest; he showed also by numerical 
examples that these waves were the same as those calculated by Stokes’s 
small-amplitude approximation; and that the value of Stokes’s parameter 
for the highest wave was approximately 0-2919. ‘This is very probably the 
actual radius of convergence. From a practical point of view, Havelock’s 
work is of no less value than Levi-Civita’s and deserves at least some 
mention. Similarly some mention should have been made of the theoretical 
work of Penney & Price (1952) on standing waves of finite amplitude, 
and of ‘Taylor’s (1953) experimental verification that the highest standing 
wave has a crest angle in the neighbourhood of 90° (as compared with 
120 for the progressive wave). 

Since the book deals only with irrotational waves, the limitations of the 
irrotational theory should be pointed out. For some years now it has been 
known that waves are not irrotational beyond the linear (small-amplitude) 
approximation. No matter how small the viscosity, some vorticity always 
spreads inward from the boundaries and results in a mass-transport 
distribution quite different from that indicated by potential theory. ‘Thus 
the irrotational theory is applicable only for a limited time after the wave 
motion is set up. 

‘The general impression left by this book is that it 1s primarily a collection 
of personal contributions, some of more interest and validity than others, 
bound loosely into one volume. ‘lhe scope of the author’s researches 1s 
wide, but the claim to have presented ** a connected account of the theory of 
wave motion in liquids, together with its applications ”’ can hardly be said 
to be substantiated when so many important topics are omitted, and others 
are treated incompletely or by passing references. In particular, the 


central position of the cnoidal wave theory in providing the link between 
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the small-amplitude approximation and the shallow-water approximation 
should have been more clearly brought out. By its very lack of unity the 
book may prove to be valuable as a collection of scattered work, although 
this is not quite what the title leads a reader to expect. 

A more serious defect is the general lack of connection between the 
mathematics and the physics. ‘This is particularly important in a field of 
applied mathematics where approximations of all kinds play an essential 
part. In at least one case—in the section on breaking waves—the author 
has been led to wrong conclusions. ‘True, there are a number of illustra- 
tions, but with the exception of chapter 11 these are not generally 
presented in such a way as to make quantitative comparisons possible. 

‘The mathematical text is well reproduced, but the photographs and 
diagrams are of poor quality in a few cases. In particular figures 6.6.1, 
10.9.1, 10.10.22 and 12.2.1 need redrawing, being rather crude sketches 
which could easily have been made more precise by use of available data. 
Many better illustrations of a solitary wave are available than figure 10.9.2. 
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Proceedings of the 5th Midwestern Conference on Fluid 
Mechanics, edited by A. M. Kurtue. Ann Arbor: ‘The University 
of Michigan Press, 1957. 388 pp. $8.00. 


This volume contains 26 of the 32 papers presented at the 5th Mid- 
western Conference, held at the University of Michigan in April 1957. 
‘The papers range over many different branches of fluid mechanics, and it 
is not possible to review the volume as a whole. Some of the papers that 
I was able to assess seemed to me to be good, although some others did not. 

Although no specific review is called for, it may be useful to notice the 
volume as another example of the growing practice of publishing the 
proceedings of conferences of all kinds. ‘There is an obvious case for 
publishing the proceedings (or some part of them) of a deliberative 
conference or symposium which has performed a job of work and has done 
something which the participants could not have done without meeting 
together. One might also argue in favour of publishing the proceedings 
of a non-deliberative conference on a restricted and important subject, 
on the grounds tat even if the papers presented could be published in the 
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regular journals it would be useful to have them all in one volume. But 
what about a conterence at which, as at that under notice, the range of 
subject matter is enormous and the organizers presumably accept for 
publication whatever papers (within certain wide limits) are offered? 

The purpose of a conference is undoubtedly to enable people to confer, 
and this it usually does well. Scientists and engineers are gregarious folk, 
as well as being internationally minded, and in these days of generous 
support from various international scientific unions and national agencies 
people are willing to go anywhere to talk about anything. Everyone has 
a good time, the equivalent of a vast amount of communication by paper 
is accomplished, and people have their outlooks enlarged by personal 
contact with a few leading spirits. ‘This is fine, and as it should be. But 
why is the legitimate business of a conference confused with the matter 
of providing a written record of individual contributions to research? 
What a man says in a talk to a conference audience is not likely to be suitable 
for publication in exactly the same form, at any rate not if it is a good talk; 
the needs of the talk and the written record are different, not only in choice 
of wording but also in structure of the account. A talk is an end in itself, 
and it would be a pity if the prospect of publication of conference proceedings 
led speakers to prepare their talks with the needs of the written record in 
mind. ‘The organizers of the Heat Transfer and Fluid Mechanics Institute 
(California) have the right idea; they issue preprints of papers to be read 
at the meetings and state that ‘‘ the preprints are made available for discussion 
only and not as contributions in tinal form to the permanent literature. 
‘The authors are tree to submit their papers for publication to any journal 
of their choice.’’ However, since these preprints are reproduced and sold 
in book form by a regular publisher, the statement does not correspond 
completely with practice. 

What puzzles me is that the practice of publishing in book form an 
almost random collection of papers read at a conference is becoming so 
common. QOn the face of it, it seems to suffer from the serious practical 
disadvantage that it does not serve the immediate interests of any of the 
parties concerned. ‘lhe publishers have on their hands a volume which is 
unlikely to sell outside libraries and is hardly an economic proposition 
without a subsidy; the editors have the thankless task of cajoling a large 
number of authors into working to a schedule and then of imposing some 
kind of uniformity on the type-scripts; and the speakers at the conference 
are obliged to prepare a complete written account of their talk, even though 
their talk may be none the better, and may even be a little the worse, for the 
existence of a straight-jacket in the form of a detailed type-script complete 
with equations; moreover, their work (if not already published elsewhere 
in slightly different form) then appears, after an unpredictable delay, in 
an irregular publication of the kind that tends to become inaccessible. 
Someone in a position of influence must think it is a good idea. Or are the 
organizers of conferences simply being too conscientious about producing 
some tangible evidence of their work ? 


G. K. BATCHELOR 
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